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Numerical calculation of flow fields about rectangular wings 
of finite thickness in supersonic flow 

Jerald Milo Vogel 

Under the supervision of Dr. E. W. Anderson 
From the Department of Aerospace Engineering and 
Dr. George Serovy from the Department of Mechanical Engineering 

Iowa State University 

The inviscid flow fields about a three-dimensional rectangular 
wing of finite thickness at angle of attack with a subsonic tip in 
a supersonic flow are determined by applying a second order finite 
difference technique to the gas dynamic equations of motion in 
their conservative form. The analysis includes a comparison of the 
second order technique with a current third order method. 

The principal objective of the study is to apply a current 
finite difference technique to the equations governing the super- 
sonic flow past a wing to obtain the variation in the gas dynamic 
variables throughout the immediate flow field. The study is 
separated into two parts. The first part deals with the comparison 
of the second order MacCormack technique and the third order Rusanov 
technique. The second part is the actual implementation of the 
numerical method to obtain the flow field about a rectangular wing 
with a 7.5 degree half-angle double-wedge cross section and a 
double-cone tip at a Mach number of 2 at 0 and 4 degrees angle of 
attack. 

The results obtained in the application of the MacCormack and 
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Rusanov techniques to the modified Burgers* equation and the gas 
dynamic equations governing the supersonic flow past a two-dimensional 
wedge indicate that the second order method MacCormack is as good 
as Rusanov's technique in terms of flow field resolution and better 
in terms of computer storage requirements and run times. 

The flow field about the rectangular wing is separated into 
three regions consisting of the forebody, the afterbody and the 
wing wake. Solutions for the forebody are obtained using conical 
flow techniques while the afterbody and the wing wake regions are 
treated as initial value problems . The numerical solutions are 
compared in the two-dimensional regions with known exact solutions. 
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NOMENCLATURE 

Coefficient matricies of the inviscid equa> 
tions of motion in the Cartesian frame 

Coefficient matricies of the inviscid equ& 
tions of motion in the ( £ , y , 0) and 
(£ , 9, 0) frames 

Speed of sound 

X-dependent conservative variables 
^-dependent conservative variables 
Constant 

Y-dependent conservative variables 
y or ©-dependent conservative variables 
Z-dependent conservative variables 
0-dependent conservative variables 
Conservative variables 
Unit matrix 

Imaginary constant (~\/- 1.0 ) 

Unit vector along streamline 

Unit vector in radial direction 

Grid indicies 

Constant 

Mach number 

Unit normal vector 

Pressure 

Differential operator 
Total velocity 
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VI 


t 

U 

U, V, w 
x, y, z 


Xj 

(f .y. *) 

(f , O, 0) 
p 


Time 

Vector of dependent variables 
Velocity components in x, y and z directions 
Independent variables in the Cartesian frame 
Angle of attack 
Shock angle 

Ratio of specific heats 

i . 

Stability parameter 
Eigenvalues 
At/Ax c 

Independent variables for body oriented 
coordinate systems 

Density 


Subscripts: 

c 

3 

k 

max 

w 

CD 


Corrected value or cone body conditions 
Mesh point location in the y, y or © directions 
Mesh point location in the z or 0 directions 
Maximum 

Conditions on the wedge 
Free stream conditions 


Superscripts: 

n 


Time or x step location 
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INTRODUCTION 

The aerodynamic analysis of the flow fields about aircraft 
capable of operating in the supersonic regime for extended periods 
of time is a formidable task. The complex geometry of such vehicles 
in conjunction with the difficulty of solving the equations govern- 
ing the aerodynamics preclude the possibility of obtaining exact 
solutions for the associated flow fields. These fundamental dif- 
ficulties have prompted the development of numerous approximate 
methods for analyzing fluid flows. One of the most common of the 
simplifying assumptions used is that the flow may be separated into 
a viscous boundary layer flow and an outer inviscid flow which ef- 
fectively determines the body pressure. This report is concerned 
with the calculation of the outer inviscid flow about a rectangular 
wing moving at supersonic speeds. 

The inviscid equations of motion governing the flow generated 
by a wing moving at supersonic speeds form a set of hyperbolic 
differential equations. Since, they are hyperbolic, the equations 
can be solved (at least conceptually) by techniques applicable to 
initial value problems. Up to the present only two such techniques 
which provide exact solutions have been applied to inviscid super- 
sonic flow problems. The first technique involves the method of 
characteristics (1). This method has been successfully applied to 
numerous supersonic flow problems. Unfortunately, the application 
of this method is a complex task due to the geometric problems 
introduced by body shape, and difficulties in determining the 
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coordinate system or systems required and the inherent way in which 
a characteristics method works. The second method involves the 
use of shock-capturing finite difference approximations of the 
equations of motion and the solution of the resulting approximate 
equations at each grid or mesh point. This attack provides a solu- 
tion for the inviscid flow throughout the flow field. The technique 
used advances the initial data through the fixed mesh, applying 
boundary conditions only at the body and in the free stream. Shock 
and expansion waves form and decay automatically without special 
treatments of any kind. On the other hand, the characteristics 
method utilizes logical numerical procedures to isolate shock waves 
and requires the application of the Rankine— Hugoniot shock relations 
across them to identify their strength and position. 

The acceptance of the shock-capturing numerical techniques is 
becoming more universal as these techniques are improved. The 
early problems associated with the precise location of the shocks 
and the tendency of the techniques to produce spurious oscillations 
in the magnitudes of the dependent variables in the neighborhood of 
the shock are gradually being overcome. Numerical calculations of 
inviscid flows based upon the full Eulerian equations have been 
carried out for a variety of supersonic problems using several 
finite difference techniques. The techniques have generally been 
first, second and, more recently, third-order. Numerous authors 
have applied the Lax ^ first-order method to fluid flow problems. 
Notable among the results obtained by these investigators are the 
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solutions for the time dependent blunt body problems obtained by 
Bohachevsky and Mates and Bohachevsky and Rubip^^ and the non- 

/ 4 \ 

equilibrium gas dynamic calculations of DeJarnette' . While the 
Lax method provides reasonable results for very small mesh sizes, 
second-order methods are being used with increasing frequency. 
Kutler^^ has recently applied a version of the second-order Lax- 
Wendroff method developed by MacCormack^^ to study flow about 
sonic-edged, conical, wing-body combinations at angle of attack. 
Results of his work show excellent agreement with conical flow 
solutions calculated using other methods and with available ex- 
perimental data. More recently, a third-order method developed 
simultaneously by Rusanov' ' and Burstein and Murin' gives im- 
proved shock and flow field resolution in certain cases. 

This study is concerned with applications of the second-order 
MacCormack^^ technique and the more recently developed third- 
order Rusanov^^ technique to some simple nonlinear problems lead- 
ing to solutions of the full Eulerian equations for flow about a 
rectangular wing moving supersonically. The material presented is 
separated into four major sections. The first section is an 
analysis of the differencing techniques under consideration as 
well as a discussion of the theoretical stability criterion based 
on amplification matrix theory. The second section is concerned 
with solutions of a one dimensional partial differential equation, 
the modified Burger’s equation. Such solutions aid in understanding 
the MacCormack and Rusanov differencing techniques and their 
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application to nonlinear hyperbolic systems. The third section is 
concerned with the application of the two techniques to a more 
realistic flow problem, the supersonic two dimensional wedge flow 
field, in am attempt to determine the technique best suited for 
the rectangular wing problem. The fourth and final section presents 
the numerical solutions for the flow fields about and in the wake 
of a rectangular, not-so-thin wing in a supersonic flow field at 
angle of attack. 
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DIFFERENCING METHODS 
Introduction 

In recent years there has been an ever increasing use of 
finite difference methods in the reduction of continuous systems 
in order to obtain solutions to complex flow problems. The dominant 
influencing factor in the development of the numerical techniques 
has been the advent of the high-speed computing machinery required 
to process data at many points in a solution field. 

To utilize a finite difference method one must first degenerate 
the continuous domain of interest to a discrete set of points 
generally termed the grid. The partial differential equations of 
motion governing the flow field are then differenced in some pre- 
scribed fashion. This results in a set of finite difference equa- 
tions which must be solved at each point in the grid subject to 
certain boundary conditions applied at the edge of the grid. 

Although the techniques appear to be elementary in nature and 
simple to apply one must be concerned with the accuracy, convergence 
and stability characteristics of the techniques. For the flow 
fields that contain shock waves one must be concerned with the 
ability of the finite difference technique to develop apparent 
discontinuities at the proper locations without producing excessive 
fluctuations in the magnitudes of the dependent variables near the 
discontinuities. In summary, improperly applied numerical methods 
may lead to extremely erroneous results. 
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Kutler^^ recently investigated a series of second-order finite 
difference techniques including Lax-Wendroff v ' methods and other 
somewhat similar methods developed by Leith ^ ^ , Fromm ^ ^ , Richt- 
meyer^^, Burstein^ 8 ^ , Strang , Gourlay and Morris ^ and 

MacCormack ^ ^ as well as the classic first-order technique by Lax ^ 
The results of the investigation indicate that in terms of ease of 
programming, storage space requirements, length of computing time 
and shock resolution and stability the method by MacCormack is 
superior and, as a result, is here to be considered for the study 
of the rectangular wing problem. 

A second technique, developed by Rusanov' ' and Burstein and 
Murin^ 8 ), will be considered to see if the more recently developed 
third order technique performs enough better so as to warrant its 
use in the rectangular wing solution. 

The following three sections present brief discussions of 
accuracy, stability and the finite difference techniques under 
consideration. 


Accuracy 

In general the errors associated with finite difference solu- 
tions may be separated into two basic types. The first type of 
error is termed truncation error and is a measure of the degree to 
which the finite difference equations actually represent the 
continuous system of equations. Truncation error may be viewed in 
terms of a set of "modified partial differential equations" which 
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is the set of equations the finite difference formulation actually 
represents (see Ref. 16, p. 38). In this investigation the "modi- 
fied equations" are studied for each differencing technique in the 
section of the report containing the technique descriptions. 

The second type of error is termed round-off error which is 
sometimes referred to as computational error. Such errors are a 
result of the discrete equations being solved exactly only up to a 
certain number of digits depending upon the particular machine used 
to obtain the solution. 

The grid point spacing effects both types of errors but in 
different fashions. That is, while decreasing the grid point spac- 
ing will generally decrease the truncation error , the resulting in- 
crease in required solution steps will tend to increase the computa- 
tional errors. Hence, one cannot always increase accuracy by de- 
creasing the mesh size. 

Stability Criteria 

One problem encountered in the application of finite difference 
schemes involves the numerical instabilities which result in error 
amplification. Unstable numerical schemes allow the growth of 
error to the extent that the true solution is "masked" yielding 
highly useless data. Hence, it is very desirable to have a means 
of predicting the parameters and their bounds which cause numerical 
schemes to result in instability. 

The stability analysis used in this paper is that outlined by 
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Richtmeyer and Morton 


(17) 


which was also utilized by Kutler 


( 5 ) 


Four restrictive conditions must be applied in utilizing the tech- 
nique. The governing equations must be locally linearized, the 
coefficients must be constant, the solution must be smooth and the 
boundary conditions must be ignored. Hence, the analysis is good 
only for regions removed from the boundaries and which are devoid 
of discontinuities. However, experience has shown that instabili- 
ties generally manifest themselves initially in the form of small 
amplitude, short wave length oscillations superimposed on a smooth 
solution in a narrow region of the solution field. Hence, the 
restrictions imposed by the stability theory may not be as prohibi- 
tive as they first appear. 

Consider the system of partial differential equations in 
conservative form given by 


E + F = 0 (1) 

t x 

where E and F are conservative variable vectors. This set of equa- 
tions can be written in the form .. 


E + A E =0 
t x 


( 2 ) 


where A is a matrix containing the Jacobian elements of F with 
respect to E. 

If the A matrix is constant one may obtain the exact solution 
by means of the Fourier series method yielding 


E = E e 
o 


-ik At ik x 
x e x 


( 3 ) 
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where E is a constant vector and k is an arbitrary constant, 
o x 

To apply the stability analysis to the numerical technique one 
introduces a row of errors along a t = constant line and observes 
the manner in which the differencing technique propagates the errors 
in t ime . The errors may be represented by a discrete Fourier series 
of the form 


F - 1 

k B k 6 


a 


(4) 


Usually, only the effects of one term (e^/^*^*) of the series is 
evaluated and a linear superposition process is utilized to evaluate 
the total error effects. If, for a fixed mesh, the error increases 
without bound as nAt o o the technique is termed stepwise unstable 
To illustrate the concept, consider the following difference 
scheme as applied to Equation (2): 


» n+1 _ =• n At A n - n 

‘ 5 7 & 2 ( j+l J-l ' 


with the boundary condition 


E ° = .i/fcA* 
D 


(5) 


( 6 ) 


The use of a separation of variables technique leads to a solution 
of the form 

i/3 3 


E ^ 11 = E x (nAt) e 


( 7 ) 


Substitution of Equation (7) into the difference equations (5) 
yields 
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E . n+1 = E . n 
3 3 


I - 


KA (e^ - 


= G E. n 
3 


( 8 ) 


where I is the identity matrix and JS is the mesh ratio. 

The matrix G in Equation (8) is termed the amplification 
matrix and the solution, given by Lomax is 



C . 
3 



(9) 


where n is an exponent and the Xj represent the eigenvalues of the 
amplification matrix. 

For the solution given by Equation (9) to remain bounded as 
n oo the eigenvalues of G must be less than or equal to unity. 
Hence, the stability criterion is given by 



( 10 ) 


To simplify the example under consideration assume that the 
set of conservative variables E contains one member. Then the 
amplification matrix reduces to 


<3 = 1 - ££ (e^ - e'^) 

In view of Equations (10) and (11) the stability requirement is 
that 


1 s + 1/ 2 A 2 s in 2 (y^Ax ) ^1 ( 12 ) 

It is noted here that the stability criterion yields stable 
values for the mesh ratio V . The results of Equation (12) indicate 
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that the difference scheme under consideration is unstable for 
any mesh ratio. 


MacCormack Technique 

MacCormack has constructed a second order predictor-corrector 
sequence for use in solving systems of partial differential equa- 
tions in conservative form. When applied to Equation (1)» Mac- 
Cormack's technique yields 

n „ n"l 
3+1 3 [ 

( 13 ) 

The tilde that appears over certain of the variables denotes the 
predicted value of that particular variable. 

To investigate the accuracy of the technique the modified 
partial differential equation is developed for a system of the form 
given by Equation ( 2 ) with A = c = constant and E = u which is the 
linear wave equation. 

u, + cu =0 (14) 

t x 

The resulting difference equations reduce to 


~ n+1 

3 


= E 


n 


At 


rJ n+1 
u . 

3 


n . . n n. 

u j - V (Uj tl - Uj ) 


u . 
3 


n+1 


- % 


fu n *v. 


n+1 , .rs n+1 

- 


r*t n+l>"l 

- Vl >J 


( 15 ) 
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The modified equation, which is the equation actually repre 
sented by Equation (13), is of the form 


u + cu = Qu ( 16 ) 

t x 


where Q is some differential operator. To evaluate Q one first 
combines the predictor and corrector equations yielding 


n+1 . .2v n 

u . a (1- 1/ )u . + 

J v D 


2 


u 


n 


3+1 


v a y +ij 

2 


u 


n 


D - 1 


(17) 


Each term is then expanded in a Taylor series about the point 
( nAt, jA?0* Partial derivatives with respect to time that are 
second order and higher are eliminated using Equation ( 16) . For 
MacCormack's technique Equation (14) reduces to 

u + cu = - cA* (l -1/2) u - C (I /-I/ 3 ) +.... 

x 6 ' ' xxx 8 ' i xxxx 

(18) 


The Qu term is represented by the right side of Equation ( 18) . 

It is noted that to second order the modified equation is 
exactly the linear wave equation as it should be since the tech- 
nique is second order. The lowest order dispersive error and dis- 
sipative error is given by the first and second terms respectively 
on the right side of Equation ( 18) . 

An interesting point to observe is that when 1/ - ± 1 the 
error terms in Equation (18) all become zero resulting in an exact 
solution. This condition is referred to as the u shift condition” 
by Kutler and Lomax who have shown that satisfying this condition 
as best possible in the nonlinear case generally yields good shock 
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capturing characteristics. 

The limits on 1/ for which the computational errors do not 
grow in an unbounded fashion may be determined by means of the 
previously discussed amplif ication matrix theory. For the linear 
wave equation the amplification matrix, which has but one element, 
is given by 


G ■ 1 - 


2 

l / cos 


(SA») - 


i j/sin 




(19) 


If the magnitude of the amplification factor, G, is not to exceed 
unity then 1 / , commonly termed the Courant number, must not be 
permitted to exceed unity. Hence, the stability bound on the mesh 
size is given by 
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In view of Equation (18) the shift condition is the maximum stable 
Courant number . 


Rusanov’s Technique 


Recent improvements in high-speed computers has resulted in 

increased interest in higher order differencing methods to improve 

flow field resolution. One of the more recent is a third-order 

(7) 

method developed simultaneously by Rusanov' ' and Burstein and 
Murin' ' . This technique, based on the Runge-Kutta method, utilizes 
a three-level predictor-corrector sequence which, when applied to 
Equation (1), is given by 
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E ( X ) - 1* (E n 

E j+h * ^ E j+i 


+ E. n ) - 1/3 
3 ' 


£ (F. n 
.x 1 j+1 


n. 


" F/ ) 




( 21 ) 





n 



) + 2F 


3-2 






- 4 


<Vi 


n 


_ n. 
+ E j.i ) 


+ 6E . 
3 


n 


+ E. n 

3-2 _ 


The last term in the third level equation is a stabilizing term 
without which the system would be unconditionally unstable for all 
values of V . 

Application of Rusanov’s technique to the linear wave equation 
given by Equation (13) yields the following modified partial dif- 
ferential equation: 


u A + cu 

t X 


-C - 4 V + l/ 3 ) U 


xxxx 


- c -Ag- (58 - 4 - 15 V 2 + 4Z/ 4 ) 


u + 

xxxxx 


( 22 ) 


It is readily apparent from the modified equation that to third 
order the linear wave equation is solved exactly as could be ex- 
pected since Rusanov's technique is third-order. The lowest order 
error term which contains the fourth derivative is dissipative in 
nature with the next higher order term being dispersive. It is also 
noted that for 8=3 and V = 1 the error terms shown on the right 
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( 19 ) 

side of Equation (21) vanish. Kutler and Lomax'’ ' have shown 
that under these conditions the shift condition is satisfied 
yielding an exact solution. 

(7) - 

In so far as stability is concerned Rusanov' has shown 
that the stability criteria are given by 



*v 2 - ^ 4 <S < 3 

which indicates that the shift condition also satisfies the 
stability requirement. It would appear that when one operates at 
mesh ratios less than unity the value of § should be set to 
most nearly satisfy the shift condition, which is a difficult re- 
quirement to meet. 
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SOLUTION OF THE MODIFIED BURGERS * EQUATION 
Introduction 

A major problem encountered in the application of finite dif- 
ference techniques is the effect of eigenvalue variation throughout 
the flow field. Stability analyses have shown that eigenvalue 
magnitudes determine stability bounds which ultimately dictate ac- 
ceptable grid mesh ratios. In addition, the modified partial dif- 
ferential equations are used to predict the best grid mesh ratios 
from an accuracy viewpoint. From previous work it has been noted 
that for best results one should operate as closely as possible to 
the upper stability bound corresponding to a Courant number of unity. 

i 

Utilizing a fixed coordinate system in which the mesh ratios are 
constant to determine a flow field in which the eigenvalues are non- 
constant precludes the possibility of operating at the best Courant 
number throughout the flow field. Hence, it is quite desirable to 
use a finite difference technique capable of good flow field resolu- 
tion through a wide range of Courant numbers. 

An investigation of the two numerical techniques under consid- 
eration is presented in this section in an attempt to evaluate 
their behavior when applied using a variety of off-design Courant 
numbers. Particular attention is given to the spreading of dis- 
continuities and oscillations of the solution near points of rapid 
change of the dependent variables. 

The hyperbolic form of the equation introduced by J. M. Burgers 
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is a valuable aid in studying the ability of a given numerical 

. - ( 20 ) 

method to produce a solution to a nonlinear equation • The 
modified Burger’s equation:: in conservative form is given by 



d* 



( 24 ) 


Kutler has successfully used this equation as an analog of the 
inviscid Euler equations and studied the solutions produced by 
various first- and second-order methods^^. Shocks and rarefactions 
which occur in the gas-dynamic solutions were simulated by introduc- 
ing discontinuities in the initial data. 

A similar procedure is followed in this section to compare 
MacCormack's technique, which Kutler found to be a superior second- 
order method, with the more recent third-order technique developed 
by Rusanov. In order to accomplish this, two discontinuities of 
different magnitudes are introduced in the initial data simulating 
shocks of different strengths with different propagation rates. 

In particular, the problem is to determine the solution of 
the modified Burger’s equation subject to the initial conditions 
shown in Figure 1 which are 



x >x 2 

x 2 >x>x 1 (25) 

X<X x 
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Since this problem represents the intersection of two discontinuities, 
the exact solution is presented in two regions. The first region 
is prior to the intersection of the discontinuities and the second 
is after the intersection (Figure 2). The exact solutions in these 

regions are 


Region 1 
u(x, t) = 0 

u(x, t) = u 2 

u(x, t) = 

Region 2 
u(x, t) = 0 

u(x, t) = ^ 




V U 2 





(26) 



Stability Analysis 

Next a stability analysis based on amplification matrix theory 
is performed on Equation (24) to determine the bounds on the mesh 
ratio At/A x for which the numerical techniques are stable. 

Equation (24) can be written in the form of Equation (2) with 
E = u and A = u yielding an expression in the form 

u t + u u x - 0 


(27) 
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which is similar to the linear wave equation with the wave speed 
equal to the quantity u. 

In view of Equations (20) and (23) the stable range of mesh 
ratios for both the MacCormack and Rusanov methods is given by 

„£<i <-> 

For those points in the solution field where At/ A x numerically 
exceeds the reciprocal of u one can expect an instability to exist. 
Numerical solutions are usually obtained by using fixed intervals 
in time and space ( At. A x )* Hence, one must search the field to 
determine the maximum of the eigenvalue, u, for this value will 
determine the largest stable At/ Ax. That is, for stable solutions 
using a fixed mesh ratio the stability criterion is given by 


c ( 29) 

^ u 

max 

It is noted here that in more complex problems the maximum eigen- 
values cannot always be determined until the numerical process is 
under way. Hence, it is not an unusual procedure to change the 
mesh ratios as the numerical process continues. 


At 

A x 


Numerical Solutions 


Two double-shock geometries were considered for each numerical 

technique. In one case u^ = 3 and u^ = 5 whereas in the other case 

u = 1 and u =5. They are termed the 5-3-0 and the 5-1-0 prob- 
2 1 

lems respectively. In both cases the discontinuities, called 
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waves, were located at x^ = 15 and x^ — 36 and were contained within 
one interval. A total of 100 intervals was used in the x direction 
with the interval size ^x equal to unity. 

The integrations in time were allowed to proceed for a total 
of 60 steps, a sufficient time interval to allow the faster moving 
wave to overtake the slower. The time interval for the integration 
was chosen to be the maximum allowable consistent with the criterion 
given in Equation (28). Under these circumstances the large ampli- 
tude wave is always being computed at the maximum Courant number 
which, according to the linear analysis, should yield the best 
possible solution. For the 5-3-0 and 5-1—0 cases the smaller 
amplitude waves are being computed at suboptimal Courant numbers 
of 0.6 and 0.2 respectively. 

The dependent variable is held constant at both ends of the 
spatial grid which provides boundary conditions for the system. 

The integration is terminated well before the waves intersect the 
boundary. 

The 5-3-0 double-shock problem is solved first. Figures 3, 

4 and 5 represent solutions using Rusanov’s technique with a 
Courant number of unity and a stability parameter, 3 , of 3, 2 and 
1 respectively. Figure 6 represents a solution using the MacCormack 
technique with a Courant number of unity. Figure 3 indicates that 
at V = 1 and S = 3 a stable solution exists throughout the field. 
The large amplitude wave is being computed at the optimum condi- 
tions while the small amplitude wave is being computed at a 
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Figure 4. Burger’s equation solution, Rusanov technique. 

5-3-0 case for = 2.0 and Courant number = 1 




Figure 5. Burger’s equation solution, Rusanov technique. 

5-3-0 case for § = 1.0 and Courant number = 1 






Figure 6. Burger's equation solution, MacCormack technique. 
5-3-0 case for Courant number = 1 
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suboptimal Courant number of 0.6. The magnitude of the stability 
parameter, § , falls within the stable range throughout the field. 

As predicted, the large wave resolution is very good with few 
oscillations occurring at the discontinuity and, at the same time, 
the discontinuity is confined to one interval. The smaller amplitude 
wave, however, is not as well behaved. The off-design Courant num- 
ber at this point in the solution field causes undesirable oscilla- 
tions to occur at the discontinuity as a result of the dispersive 
terms in the modified equation. The spreading of the discontinuity 
caused by the dissipative terms in the modified equation is not 
extremely significant at this Courant number since the discontinuity 
remains captured in essentially one interval. Figure 4 shows the 
results for 1 / - 1 and § ** 2. For the large wave § is in the un- 
stable range according to the linear theory. Although it appears 
that the actual numerical solution is stable an excessive number of 
large amplitude oscillations occur yielding a highly undesirable 
solution. For the lower amplitude wave the § value is in the stable 
range. The oscillations at the discontinuity are fewer than the 
previous solution indicating that perhaps the stability parameter 
should be set at less than the maximum value for the best results 
at off-design Courant numbers. In both cases the discontinuities 
remain isolated between essentially two grid points. Figure 5 in- 
dicates a solution for which 1 / = 1 and § = 1. For the large 
amplitude wave the value of S is far outside the stability range 
and, as a result, an instability occurs at this point in the solution 
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field. For the small amplitude wave the § parameter is slightly 
outside the linear stability range yielding a solution with excessive 
oscillations but apparently marginally stable. Figure 6 shows the 
MacCormack solution at a Courant number of 1.0 for the large wave 
and a Courant number of 0.6 for the small wave. Essentially no 
oscillations occur at the large wave while even at the off -design 
Courant nu mb er very few oscillations occur at the small wave. Both 
discontinuities remain isolated between two grid points. 

Figures 7, 8 and 9 represent a Rusanov solution to the 5-1-0 
problem with IS = 1 and for values of 8 of 3, 2 and 1 respectively. 
The same general trends occur with the addition of substantial 
amounts of discontinuity spreading at the lower wave which is 
being computed at a Courant number of 0.2. Figure 10 shows the 
MacCormack solution which contains, overall, fewer oscillation 
problems and a lesser amount of discontinuity spreading. 

On the basis of the information obtained from the solution of 
Burger’s equation it appears that it is desirable to use the Mac- 
Cormack technique rather than the Rusanov technique. Shock resolu- 
tion and over- and under-shoot characteristics are better over the 
range of eigenvalues considered. At the lower Courant numbers the 
shock spreading with the MacCormack technique appears to be less 
severe than that resulting from the use of the Rusanov technique. 

In addition, the computer storage and computation time requirements 
are significantly lower for the MacCormack technique. 



Figure 7. Burger's equation solution, Rusanov technique. 

■5-1-0 case for q = 3.0 and Courant number = 1 



Figure 8. Burger’s equation solution, Rusanov technique. 

5-1-0 case for § = 2.0 and Courant number = 1 
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NUMBER OF TIME STEPS 


Jl 

Figure 10. Burger’s equation solution, MacCormack technique. 
5-1-0 case for Courant number = 1 
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TWO-DIMENSIONAL WEDGE FLOW FIELDS 
Introduction 

To further develop the comparison between the numerical tech- 
niques under consideration a study of the two-dimensional wedge in 
supersonic flow is undertaken. The equations of motion governing 
the flow about a two-dimensional wedge with the same half-angle as 
is encountered at the leading edge of the three-dimens ioftal rectan- 
gular wing are solved at a Mach number of 2 using both numerical 
techniques. Then a comparison with available exact solutions is 
made to evaluate the performance of the techniques. Of particular 
interest is the capability of the techniques to develop crisp shocks 
in the proper locations as well as minimize the number of oscilla- 
tions of the dependent variables in the neighborhood of the shock. 

Three different approaches may be used to obtain a numerical 
solution of the wedge equations of motion. 

In the first approach the complete unsteady equations of motion 
are integrated subject to boundary conditions dictated by the wedge 
geometry. Since the equations are hyperbolic in the time variable, 
the problem generated is of the initial value type. The integration 
of these equations in time proceeds from an appropriate set of 
initial data and is terminated once the flow variables reach a 
steady state condition. 

The second approach is one used by Kutler^. The full-blown 
equations of motion are reduced to a set of steady equations by 
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setting all time derivatives equal to zero which would be the situa- 
tion in the case of the steady state solution. The resulting equa- 
tions are hyperbolic with respect to the x-coordinate (the coordinate 
most nearly aligned with the. flow direction) as long as the x-compo- 
nent of velocity remains locally supersonic. Again, the system 
reduces to an initial value problem and can be integrated in the x- 
direction starting from an appropriate set of initial data and sub- 
ject to boundary conditions dictated by the wedge geometry. Since 
the flow is conical in nature the flow variables in the solution 
are constant along rays from the origin, a condition which serves 
as the convergence criterion in the numerical process. It is readily 
apparent that the x-coordinate in the steady equations is quite 
analagous to the time coordinate in the set of time dependent 
equations. 

The third approach is one that has been used by Anderson and 

( 21 ) 

Vogel' . The full-blown equations of motion are transformed 
from a Cartesian coordinate system to a polar coordinate system in 
which one of the coordinates (r) is the distance along a ray from 
the origin. Again, the set of equations is hyperbolic in time. 

Since the flow is conical, the steady state values of the r deriva- 
tives are initially set equal to zero yielding a simple set of 
hyperbolic equations containing one less independent variable than 
in the case of the full-blown set. The equations are integrated 
in time starting from an appropriate set of initial data. The 
solution is realized when the dependent variables no longer change 
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with time. 

In general, approaches 2 and 3 require less computer storage 
and computation time them is the case in approach 1, since the 
equations of motion contain one less independent variable than 
those of approach 1. However, approach 1 is probably more versatile 
since it is not dependent upon the conical flow requirement. 

Approach 2 is adopted for the investigation undertaken in this 
paper. The numerical solutions for the two-dimensional wedge as 
well as the leading edge of the rectangular wing are obtained 
through numerical integration of the steady equations. 

Steady Equations of Motion 

The two-dimensional wedge flow equations of motion in a 
Cartesian body-oriented coordinate system for a steady, inviscid, 
nonheat-conducting and adiabatic flow are given by^ 5 ^. 

d{P u) + d(/° v ) = o 
d * 3 y 

P +Pu 2 ) + d(Puv) _ Q 

a x dy 

(30) 

$(P UV ) + d( p +Pv 2 ) _ o 

d x 3 y 

p = f 1 + ^ 2 + v2 j] 

These equations are the continuity equation, x-, and y-direction 
momentum equations and the integrated form of the energy equation 
which is usually referred to as Bernoulli’s equation. The dependent 
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variables (yO» u, v, p) in the equations are in a dimensionless 
form. The nondimens ionalizing parameters for the pressure, density 
and the velocity components are gamma times the free stream stagna- 
tion pressure, the stagnation density and the stagnation speed of 
sound respectively. 

Three partial differential equations of motion are in the 
conservative form 


-£l + j2£ = 

d x By 


o 


( 31 ) 


where E is a vector whose elements are conservative variables 
given by 



and F is a vector whose elements are conservative variables 
given by 
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Exact Solution 


The exact solution to the two-dimensional supersonic wedge 

( 22 ) 

is clearly presented in the text written by Liepmann and Roshko' 

The flow field over the wedge surface behind the shock is uni- 
form and is in the direction of the wedge surface. The wedge 
surface pressure is given by 


P - P a 
w 

^co 




(M co 2s±T)2 /3 “ 1 ) 


( 34 ) 


and the shock angle can be determined using the equation 

2 . 2 / 


tan 0 = 2 co 


a 


M m sin 




M, 


QO 


\y*cos 2 /3 )+2 


( 35 ) 


Two values of shock angle (jQ) satisfy the equation. The smallest 
value, the proper solution, represents the weak solution for at- 
tached shocks. 


Wedge Coordinate System and Resulting Grid 


The Cartesian coordinate system used for the wedge flow field 
analysis is aligned such that the x-axis is in the direction of the 
wedge surface with the y-axis normal to the surface. Hence, the 
wedge surface is a constant-coordinate surface, a highly desirable 
situation with regard to boundary conditions. Application of 
boundary conditions for body surfaces in nonaligned coordinate 
systems can be extremely difficult and, at times, may present 
stability problems. 
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The grid system for the numerical techniques generated by 
the coordinate system is shown in Figure 11. Since only the upper 
surface flow field is to be considered the lower wedge surface can 
be ignored. The existence of the sublayer is a result of the ap- 
plication of boundary conditions and is discussed later. It is 
also noted that for this coordinate system the free stream velocity 
vector is canted with respect to the x-axis. 

A second coordinate system used in the two-dimensional wedge 
analysis but not presented in this section is depicted in Figure 
12. The system is termed a "semi-polar 1 ' system. As with the 
Cartesian system, the wedge surface is a constant -coordinate surface. 

Numerical Solution Technique 

The integration in the Cartesian system is initiated using free 
stream values of the flow variables at each grid point in the y- 
direction along the x = 1 line. This is commonly termed an 
impulsive start. The integration in the x-direction continues 
until the x = 2 line is reached. At this point in the numerical 
process Kutler's stepback procedure is impelemented^ 5 ^ . This allows 
the integration to be re-started at x = 1 with updated initial data 
taken from the x = 2 line. In Figure 13 it is noted that grid 
points numbered 2, 4, 6, . ..., m along the x = 2 line are on the 
same rays from the origin as the grid points 2, 3, 4, ... 
m/2 + 1 along the x = 1 line. Since the flow is conical, 
data may be shifted along the rays from the x = 2 to the 




Figure 11. Coordinate system and grid system for 2-D wedge 









41 


x = 1 grid points generating new initial data. For the scheme to 
work properly the grid points m/2, m must be outside the 

shock wave in the free stream. Solution convergence occurs when 
the initial data generated by successive stepbacks becomes constant. 

It is noted that for the semi-polar system depicted in Figure 
12 the stepback procedure becomes somewhat meaningless since the 
corresponding grid points at all x = constant stations lie on rays 
from the origin. Hence, the integration can proceed in the x-direc- 
tion until no change occurs in the dependent variables along the 
rays . 

Boundary conditions must be specified at the outer grid points 
and at the sublayer grid points since these are not integrated 
points. The dependent variables at the outer grid points are more 
easily handled. The grid is set up so that the outer edge is always 
in the free stream resulting in known constant boundary data. The 
grid points along the sublayer present a more difficult problem. 

The normal velocity component at the surface of the wedge must 
vanish since flow cannot pass through the surface. To satisfy this 
condition the normal velocity component is treated as an odd func- 
tion at the body surface. That is, the sublayer value of the normal 
velocity component is set equal to the negative of the normal velo- 
city component one layer above the surface. The values of the re- 
maining dependent variables along the sublayer are evaluated using 

( 3 ) 

the reflection technique as used by Bohachevsky and Rubin' ' . The 
basic assumption used in this technique is that the normal derivatives 
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of the dependent variables vanish on the surface. Hence, the 
variables are treated as even functions on the wedge. That is, 
the dependent variables along the sublayer are set equal to their 
values one mesh point above the body. In view of the exact solu- 
tion the reflection technique is exact for wedge flow in that the 
normal derivatives are zero. 

Numerical Solutions 

Equations (30) were integrated using the Iowa State University 
IBM 360-65 computer system for both Rusanov's and MacCormack's 
methods for a wedge with a 7.5° half-angle at a Mach number of 2. 
Two mesh ratios ( Ax/Ay) were used. One at 1.272 which is near 
the experimental maximum for stability as determined by Kutler 
and the other at l.O^ 5 ^. Three values of the stability parameter, 

§ , associated with the Rusanov technique at a mesh ratio of 1.0 
were used to assess the effect of S on the solution. In all cases 
the grid points in the y-direction consisted of a total of 30 mesh 
points. 

Figure 14 shows the pressure distribution normal to the wedge 
surface using Rusanov's method at the lower mesh ratio of 1.0 for 
S values of 1.0, 2.0 and 3.0. The solution for a 3 of 1.0 is 
distinctly inferior to those obtained for 3 values of 2.0 and 3.0. 
This is a result of the excessive overshoots and undershoots in the 
vicinity of the shock. The solution for 3 = 2.0 appears to yield 
a crisper shock than the solution for 3 = 3.0 as well as lower 
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Figure 14. Pressure distribution in a direction normal to 
wedge surface. Rusanov technique for a two- 
dimensional wedge with half-angle 7.5° at a Mach 
number of 2 for Ax/Av =1*0. x = 1.0 
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amplitude oscillations in the free stream after the shock is 
encountered. 

Figure 15 shows a similar pressure distribution using the 
MacCormack technique for Ax/Av = 1.0. Although there is some 
overshoot from the shock layer side, the behavior for the free 
stream side is very good with no oscillations occurring. 

Figures 16 and 17 show the pressure distributions for the 
Rusanov and MacCormack methods, respectively, at the higher mesh 
ratio near the experimental stability bound. The Rusanov stability 
parameter, § , was set equal to 3.0 which, according to linear 
stability theory, is the only stable value when the maximum mesh 
ratio is used. The MacCormack technique develops a crisp shock 
with few oscillations as the shock is encountered on either side. 

The Rusanov solution, on the other hand, develops a fairly crisp 
shock but exhibits excessive oscillations on the free stream side 
of the shock. While decreasing the value of S improves the flow 
f ield behavior in the free stream, the shock layer portion of the 
solution becomes less well behaved in this case. 

In all cases, the shock is properly located and the magnitudes 
of the dependent variables are correct. 

There is about a 30 per cent savings in computation time as 
well as a substantial decrease in storage requirements using the 
MacCormack technique. Based on these criteria as well as the re- 
sults discussed above, it would appear that the MacCormack technique 
is superior for solving the two-dimensional wedge flow problem in 
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Figure 15. Pressure distribution in a direction normal to 
wedge surface. MacCormack technique for a two- 
dimensional wedge with half-angle 7. 5° at a Mach 
number of 2 for Z\x/ Ay =1.0. x = 1.0 
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Figure 16. Pressure distribution in a direction normal to 
wedge surface. Rusanov technique for a two- 
dimensional wedge with half-angle 7.5° at a Mach 
number of 2 for Ax/ Av = 1.272 and ^ = 3.0. 
x 
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Figure 17. Pressure distribution in a direction normal to 
wedge surface. MacCormack technique for a two- 
dimensional wedge with half-angle 7.5° at a Mach 
number of 2 for Ax/ Ay = 1.272. x = 1.0 
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the range of Mach numbers and wedge angles examined here. 

Additional experiments were performed using the semi-polar 
coordinate system described earlier. Although these solutions 
are not discussed here, they proved to be satisfactory. The major 
difference noted was that computation times required to reach a 
solution were longer. 
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THE RECTANGULAR CONE-TIPPED WING 
Introduction 

The primary objective of this investigation is to develop, 

by means of finite difference techniques, the flow field about a 

( 23 ) 

body in supersonic flight. Kutler and Lomax' have already 
investigated a variety of wings, bodies and their combinations 
including two-dimensional wedges, two-dimensional and axisymmetric- 
nonconical bodies, cones, planar delta wings and delta wings with 
dihedral mounted on. conical bodies. The latter two studies were 
restricted to supersonic leading edges. His results obtained 
using shock capturing finite difference techniques agree well with 
the method of characteristics solutions and available experimental 
data. 

The body considered in this paper is a three-dimensional 
rectangular wing with a symmetric double wedge cross section to 
which is attached a double cone tip as shown in Figure 18. The 
cone half-angle is chosen to be less than that of the tip Mach 
cone yielding a subsonic tip. Hence, the upper and lower surfaces 
are not independent as they are in the case of supersonic edges. 

The flow field about the body can be separated into three 
distinct regions as indicated in Figure 19. 

Region I contains the forebody flow field which begins at the 
cone-wedge vertex and ends at the mid-chord point. The wing aspect 
ratio and free stream Mach number are chosen such that the Mach 
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cones from the wing tips do not intersect on the forebody. Hence, 
the flow in this region is entirely conical. That is, flow variables 
along rays from the vertex of the cone are constant. The flow field 
solution is generated using conical flow methods similar to those 
discussed in the section on 2-D wedges. 

Region II contains the afterbody which begins at the mid-chord 
point and ends at the aft cone-wedge vertex. The partial differen- 
tial equations governing the flow in this region represent an 
initial value problem that is solved using initial data generated 
in the forebody solution. That is, the flow variable magnitudes in 
the plane containing the base of the forebody are the initial data 
for the afterbody. 

Region III contains the wake behind the wing which begins at 
the aft cone-wedge vertex and extends downstream indefinitely. 

Again, the problem in this region reduces to an initial value 
problem with initial data generated using the afterbody solution. 

The sections that follow contain discussions of the problems 
and the solutions associated with each of the three regions. Dis- 
cussion topics include equations of motion, coordinate systems, 
stability, boundary conditions, solution techniques and results. 

Equations of Motion 

The basic flow equations that govern a supersonic flow are 
given by the conservation of mass, energy and momentum. The first 
two equations are scalar while the last is a three-component vector 
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equation yielding a total of five scalar equations of motion. For 
a steady, inviscid, nonheat-conducting and adiabatic flow these 
equations in vector notation are given by, respectively 

V* (/Oq) = 0 

q • VH t = 0 (36) 

V(q 2 /2) + (V x ^) X S + VP'/O = o 

For a Cartesian coordinate system Equations (36) may be written in 
the scalar form as follows: 

Conservation of mass 

</0«) x + (,Gv) y + (p w) z = 0 

x-direction momentum 
(p + /Du 2 ) x + (/Ouv) y + (yOuw) z = 0 

y-direction momentum (37) 

(/° UV >x + < p + = 0 

z-direction momentum 
(/0 UW ) X + (/Ovtfy + (P + /Dw 2 ) z = 0 

Energy equation 

/° r V-l , 2 . 2 2.1 

P = y - [1 - ■’ 2 — (u + v + w )J 

The dependent variables (p, p , u, v, w) in Equations (37) are in 
a dimensionless form. The nondimensionalizing parameters for the 
pressure, density and velocity components are gamma times the free 
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stream stagnation pressure, the stagnation density and the stagna- 
tion speed of sound, respectively. 

It is noted that only four of Equations (37) are partial dif- 
ferential equations. The fifth equation, the energy equation, is 
used in its integrated form to simplify the integration procedure. 

Equations (37) said equivalent equations in other coordinate 
systems can be written in the general form 

E + F + G + H = 0 (38) 

X y z 


where E, F and G are vectors whose elements are conservative 
variables given by 


E = 


/° u 

n. r 

p + p* 

yOuv 

yOuw 


\ 


(39) 


F = 



(40) 


G = 



p + pw / 


(41) 
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The vector H represents the non-homogeneous portion of Equation (38) 
and is identically zero for the Cartesian form of the equation of 
motion. However, it is nonzero in the equations of motion associated 
with the forebody and afterbody flow because the Cartesian system 
is not employed in these regions, only in the wake region. , The 
coordinate systems and associated equations of motion for regions 
other than the wake are developed and discussed in subsequent sections. 

Numerical Technique 

Based on the analyses in previous sections concerning the 
modified Burger’s equation and the two-dimensional wedge flow solu- 
tions, the MacCormack technique was chosen over the Rusanov tech- 
nique. The criteria used to make the comparison were solution time 
and computer storage requirements as well as flow field resolution. 

Over the range of eigenvalues considered the MacCormack technique 
produced crisp shocks with few oscillations on either side of the 
shock, generally better than the Rusanov technique. Solution times 
were, on the average, thirty per cent less with a substantial de- 
crease in the computer storage requirement. 

The conclusions reached thus far may be somewhat misleading in 
that there appear to be situations in which the Rusanov technique 
is quite superior to that of MacCormack. Anderson and Vogel 
have investigated the shock reflection problem in which a shock wave 
intersects a flat surface resulting in a second reflected shock. 

The equations of motion governing the flow are the normal fluid flow 
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equations discussed in this paper. One encounters a wide variation 
in eigenvalue magnitudes throughout the flow field in this situation. 
Hence, utilizing a constant mesh ratio in applying the numerical 
technique causes a portion of the solution to be developed at a very 
low effective Courant number. In those cases where the eigenvalues 
vary as much as a factor of ten, Rusanov’s technique produces 
distinctly superior solutions. However, no reflected shocks are 
encountered in the rectangular wing problem and, as a result, such 
severe eigenvalue variations do not occur throughout the flow field. 

The MacCormack predictor-corrector equations as applied to 

I 

Equations (38) are given by 


/V/ 

E 


j.“ 1 -“ J . k n -^<w n " d . k n » 


y ,k+i 


n 


n+1 




= \ 


[ _ n 

E . , + E . . 

J »k J,1 


n+1 


G. ") - H. ” Ax 

J »k ' 3,k ^ 


n+1 * n+1. 

(F j,k " F j-l,k > 


(42) 


x ,C* / n+1 

S (G j,k 


G, , n+1 ) - H, , n+1 Ax 


D »k-l 


j »k 


The tilde that appears over certain of the variables denotes the 
predicted value of that particular variable whereas n, j and k are 
the indicies associated with the x, y and z directions respectively 
and serve to define the location of the grid points throughout the 


flow field 


/ 



Stability Considerations 


The quality of solutions obtained using Equations (42) depend 
to a great extent upon the magnitudes of the mesh ratios /\x/ /\y 
and Z\x/ . Operation at mesh ratios outside the stable range 
leads to divergence whereas values well below the maximum stable 
values lead to poor flow field resolution in the neighborhood of 
the shock. Hence, it is quite desirable to have a priori knowledge 
of stable ranges of the mesh ratios in setting up the flow field 
grid. 

Kutler and Lomax present criteria based on amplification 
matrix theory to theoretically predict stability bounds in multi- 
dimensional problems. To utilize the analysis one must know the 
eigenvalues of the coefficient matricies of the gas-dynamic equa- 
tions of motion. In his work Kutler^ ^ developed the coefficient 
matricies for the equations of motion in a Cartesian system and 
determined the associated eigenvalues. For convenience, his work 
is outlined in Appendix A. The stable ranges for the mesh ratios 
are given by 



( 43 ) 


where 


and 


max 


B 


represent the maximum eigenvalues of the 


max 
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A said B matricies respectively. One can obtain approximate values 
for the maximum eigenvalues to be used in conjunction with Equations 
(43) to set up a grid in which the mesh ratios are near their maximum 
value yielding near-optimum numerical results. 

Region I Flow Field Analysis 

Forebody geometry and coordinate systems 

That portion of the body contained between the leading edge 
of the wing and the shoulder at mid-chord constitutes the forebody. 
The geometry of the forebody and the associated Cartesian coordinate 
system are depicted in Figure 20. As noted in Figure 20, the fore- 
body can be separated into two parts. The first part, the wing 
proper, consists of a wedge with half-angle 0^ while the second 
part, the wing tip, consists of a half-cone having the same half- 
angle © c as the wedge. Hence, a smooth transition is made from 
the wedge to the cone with no discontinuities in surface slope. 

Since the wing is not cambered, the forebody cross section is 
symmetric with respect to the chord plane. 

The Cartesian coordinate system associated with the forebody 
has the origin at the apex of the cone. The y-axis is perpendicular 
to the plane of symmetry, the x-axis is in the chord plane (the 
plane of symmetry) and the z-axis extends along the span of the 
wing. The positive directions are as shown in Figure 20. 

It has been noted earlier that a coordinate system in which 
the body can be described by a constant coordinate surface is highly 
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desirable. For the forebody this is most easily accomplished using 
two systems, one for the wing proper and one for the half-cone. 

The coordinate system used to describe the half-cone and the 
corresponding flow field region is shown in Figure 21. The coordinate 
0 defines the angle between the interface plane (the plane containing 
the intersection of the wedge and cone) and the cone meridian plane 
containing the point ( £ , 0, 0 ). The coordinate 0 defines the angle 
between the cone axis and the radius to the point of interest. The 
coordinate £ is simply the x position of the point. 

The coordinate system associated with the forebody wing proper 
and corresponding flow field is shown in Figure 22. Only the upper 
half of the body is shown. The coordinate 0 represents the angle 
between the interface plane and the plane normal to the chord plane 
and containing the radius to the point (£, ®* 0 ) • The coordinate 
0 represents the angle between the chord plane and the projection 
of the radius on the interface plane. The coordinate £ is the x 
position of the point. 

The grid system generated by the coordinate system is defined 
by the intersection of a set of 0 = constant planes, a set of 0 = 
constant surfaces, one of which defines the body, and a £ = constant 
plane. The grid in a typical £ = constant plane is depicted in 
Figure 23. The grid points associated with the wedge in the inter- 
face plane are identical with those of the half-cone. It is noted 
here that corresponding grid points in successive ^ = constant 
planes lie along rays from the common origin of the coordinate 
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Figure 22. Coordinate system used to describe the wing proper of the 
forebody and corresponding flow field region 
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systems. Hence, the stepback procedure becomes somewhat trivial 
in the forebody solutions. 


Forebodv flow equations of motion 

To obtain flow field solutions in Region I the flow equations 
must be transformed from their Cartesian form given in Equations (37) 
to the new systems described in Figures 21 and 22. The mechanics 
of the transformations are presented in Appendix B. 

The transformed equations of motion for both the wedge proper 
and the half-cone for the forebody region as given by Equations (B4) 
and (B5) are of the general conservative form 


E^ ' + F q ' + G^' + H* = 0 (44) 

_i _ t —i 

The conservative variable vectors E , F , and G as well as the 
nonhomogeneous term H for the wedge proper are given by 


e' = f E 

—i . — 2 — 

F = -sinGcosQE + cos ©F 

G = -sin^cos^E + cos 

H = (-2s in Q - sin 0 + cos 0)E 

+ 2sin©cos©F + 2sin0cos£&3 

and for the half-cone are given by 

E = £ tan©E 


(45) 

(46) 

(47) 


(48) 


(49) 
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I MB O — MM 

F = -Esin 0 + Fsin©cos©cos0 + Gsin©cos©sin0 
G = -Fsin0 + Gcos0 

H » (2sin©cos© - tan©)E 

— o — 2 

+2Fcos0sin © + 2Gsin ©sin0 

where the vectors E, F and G are the conservative variables 
associated with the standard Cartesian form of the equations of 
motion and are defined in Equations (39), (40) and (41). 

The scalar components of Equations (44) represent the conserva- 
tion of mass and the x, y and z direction momentum equations. The 
energy equation is not included in this set since it is used in the 
integrated form as given by the last of Equations (37). 

Evaluation of gas dynamic variables from conservative variables 
In the integration process, the set of predictor-corrector 
equations developed by MacCormack are used to determine only the 
numerical values associated with the scalar components of the vector 

t t 

E . Each time that E is updated along the integration path in the 

t t 

-direction the remaining conservative vectors F and G as well 

— i 

as the nonhomogeneous term H must be numerically evaluated. Since 

—i i i 

F , G and H are functions of p, , u, v, w one must extricate 

—i 

from E the gas dynamic variables. 

In the coordinate systems used to describe the various regions 
of the wing and surrounding flow field the variables E are simple 
functions of only the Cartesian counterpart E. Hence, it is 



(50) 

(51) 

(52) 




65 


convenient to evaluate the gas dynamic variables using a two-step 
process. In the first step the scalar components of E are evaluated. 
Secondly, the gas dynamic variables are evaluated using the inverse 
of the relationships E (p,^0 » u, v, w) and the integrated form of 
the energy equation. The explicit forms for the variables p, ^ , 
u, v and w are developed in the following work. 

In view of Equation (39) the scalar components of E are given 


by 

E 1 = yO u ( 53 ) 

E 2 = p + yOu 2 ( 54 ) 

E 3 = yOuv ( 55 ) 

E^ = yQUW (56) 

Dividing Equations (55) and (56) by (53) yields 

v = E 3 /E 1 (57) 

w = E^/E 1 (58) 


Combining Equations (53), (54) and the last of Equations (37), 
the energy equation, yields 



(7+ i) 
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The positive sign is used since flow is supersonic throughout the 
flow field. Equation (53) may then be used to evaluate density as 
follows, 

P * E i/“ 

and the energy equation for pressure 

P = ~y~ (! - q 2 ) 

The relationships between E and E for the wedge proper and 
the half-cone in the forebody region are given by Equations (45) 
and (49) respectively. 

Initial and boundary conditions 

In the application of MacCormack • s technique the grid points 
existing on the boundaries of the grid system depicted in Figure 
23 are not integrated points. Hence, a set of boundary conditions 
must be developed to specify the values of the gas dynamic variables 
along the upper and lower two-dimensional wedge boundaries, the 
outer free stream boundary and the upper and lower wing surfaces. 

The conditions along the outer free stream boundary are most 
easily specified since the grid is always made large enough that 
the outermost grid points always lie outside the shock in the free 
stream. Thus, the gas dynamic variables at these locations retain 
their constant free stream values. 

The upper and lower two-dimensional wedge boundaries are 
placed toward the center of the wing well outside the Mach cones 


(60) 


(61) 
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emanating from the half-cone ver tides of the wing tips. Hence, 

these grid points exist in the two-dimensional wedge regions of the 

wing flow field for which exact solutions are known. The gas 

dynamic variables along these boundaries are frozen at the two- 

dimensional wedge flow values dictated by the exact solution. The 

dimensionless shock layer pressure is given by Equation (34), 

( 22 ) 

the dimensionless density is given by v ' 

A = <y *1) 2 sin^ /3 ’ 

£ ~ (7-1) Moo 2 Sin 2 /? ' + 2 

and the dimensionless rectangular velocity components are given 
by 


u w (y- 1 ) M co 2 sin 2 /? ' 4 - 2 


(7+1) M a, 2 sin /? 

v u 

= — w— tan © 

q « <5® w 

w =0 
w 


(63) 

(64) 

(65) 


where /? ' = ^ upper + CC for the upper surface and /J = /? lower - CC 
for the lower surface (see Figure 24). Equation (35) defines the 
shock wave location for the upper and lower surfaces in the two- 
dimensional regions with © = Q w - CC and © » © + 'C C for the upper 

and lower surfaces respectively. 

The boundary conditions that must be specified along the 
surface of the wing in each A = constant plane are somewhat more 
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complex. Although a variety of techniques has been developed many 
of the procedures are fairly difficult to implement and sometimes 
quite costly in terms of computing times. Also some yield question- 
able results. 

Abbett ^ 24 ^ has reviewed and compared many of the procedures 
currently in use for calculation of surface boundary points. In- 
cluded in the survey are reflection techniques, explicit and implicit 
differencing using one-sided derivative approximations, characteris- 
tics techniques and techniques utilizing extrapolation from interior 

f 24 ^ 

points to the boundary. In addition, Abbett' ' has developed a 
new scheme in conjunction with MacCormack’s differencing technique 
to evaluate the gas dynamic variables along the wall. The method 
is analytically simple, easy to incorporate, computationally fast 
and satisfies an entropy condition on the body surface that the 
other techniques do not. Although the basic Abbett technique does 
not appear to give good results for bodies having high curvature, 
a slight modification yields a very usable scheme which is applied 
in this study. The following paragraphs include an explanation of 
the basic Abbett technique as well as the required modifications. 

The basic Abbett boundary condition scheme is a two-step 
predictor-corrector sequence in which the prediction step consists 
of the original MacCormack predictor. All gas dynamic variables 
on the body are evaluated in this step. In general, the predicted 
velocity at the body will not be parallel to the surface. The 
corrector step, then, consists of an application of a simple 
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isentropic flow expansion or compression, whichever is necessary, 
to turn the flow parallel to the surface. Then the body pressure 
is corrected by the amount dictated by the expansion or compression 
angle and the body density is determined using the surface entropy 
value as well as the corrected pressure. The modulus of the velo- 
city on the body surface is corrected by using the energy equation. 

The angle ( 3 ) through which the flow must be turned after the 
predictor step to align the flow with the surface is given by 


8 . sin" 1 Li! 

% 


( 66 ) 


where the subscript p denotes predicted value and ft is the Unit 
vector normal to the surface. 

Now that the turning angle for the expansion or compression 

has been evaluated, the corrected pressure may be determined using 

( 25) 

a truncated form of the Prandtl-Meyer function' ' given by 


^c 

P 

P 


i.ms + y 

Vm -l 


M 


(V-H)M 4 -4(M 2 -1) 

4(M 2 -!) 2 


(67) 


The subscript c denotes corrected value and the parameter M denotes 
the predicted value for surface Mach number. It is noted here that 
for positive § quantities the flow must be expanded and for negative 
values the flow must be compressed. 

Since the surface entropy is known and constant the body 
density may be corrected using 
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P c = P (P c / P)^ < 68 > 

where the barred quantities are known constants. 

The corrected velocity modulus from the energy equation is 
given by 


( 


y p c n 

-> t 1 “ ^ 




(69) 


whereas the proper velocity direction is defined by means of the 
unit vector 


q \% “ ( s * "I 


(70) 


Equation (70) is developed by removing the normal component of the 
predicted velocity and normalizing the resulting vector which 
effectively defines a unit vector parallel to the surface. The 
combination of Equations (69) and (70) yields the corrected 
velocity 



(71) 


from which the rectangular components u c , v c and w c may be deter- 
mined. 

Ferri^^ has shown that for conical flow the entropy along a 
streamline remains constant. Hence, the body entropy can be defined 
once the body streamline is identified. Since the streamlines that 
wet the body surface emanate from the two-dimensional regions outside 
the tip Mach cone, the entropy values can be obtained from the exact 



72 


2-D wedge solutions. For the single -of -at tack case the upper and 
lower surface shock waves in the 2-D regions are of different 
strengths yielding different body entropies for the upper and lower 
surfaces. Ferri^^ also shows that vortical singularities can, in 
fact, occur at those points in the flow field where both the normal 
velocity component (normal to the radius from the origin) and the 
crossflow velocity component vanish. Since the normal component 
of velocity is zero everywhere on the body surface a vortical 
singularity on the body can occur only where crossflow stagnates. 
Hence, to apply Equation (68) one simply uses upper surface 2-D 
values of p and yO along the upper surface 3-D region until the 
crossflow stagnation point is reached. Then the lower surface 2-D 
values must be used. 

For bodies having high curvature the prediction step in the 
basic Abbett technique appears to misalign the flow on the body 
surface. The Abbett corrector then continuously compresses (shock 
layer region) or expands (expansion region) the flow throughout the 
integration process resulting in very large or near zero body 
pressures respectively. To remedy the situation the reflection 
technique is used after the prediction step to produce more realis- 
tic body flow variables for the Abbett correction step. 

A sublayer is added to the grid system. This yields a set 
of mesh points below the body surface as depicted in Figure 23. 

The sublayer values of the flow variables are estimated by treating 
the normal velocity component as an odd function and the radial and 
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crossflow components as well as the pressure and density as even 
functions at the body surface. The normal, radial and crossflow 
scalar velocity components are given by 

q« = q * n 

q r = i • * r ( 72 ) 

q c = i - % “ ^r 


respectively, where i is the unit vector in the radial direction. 
The resulting equations defining the sublayer flow variables are 
given by 


\ = -% 


(73) 


where subscripts 1 and 3 denote sublayer and superlayer values 
respectively. 

Once the sublayer values are known the body surface grid 
points may then be evaluated using the MacCormack corrector. The 
resulting body flow variables are used as inputs to the Abbett 
boundary condition routine which satisfies the body entropy condition. 
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The integration procedure is started impulsively. That is, 
the gas dynamic variables at the interior grid points are initially 
set at free stream values. Although different initial data may be 
used to initiate the integration the impulsive start seems to be 
quite convenient. As the integration proceeds from one £ plane to 
the next, the shock wave moves from the body surface out into the 
flow field to its proper location. 


Choice of grid system 

In view of the stability and accuracy considerations presented 
earlier the forebody grid should be chosen such that the mesh ratios 
a£/a© and a£/a a are near their upper bounds for stability on 
both the wedge proper and the tip half-cone. Failure to do so may 
yield poor solutions in portions of the flow field. Of greatest im- 
portance is the a£/a © ratio since the shock wave is encountered 
in the © direction. In previous work dealing with the 2-D wedge 
solutions it was noted that the most disastrous effects of suboptimal 
Courant number operation occur in the vicinity of the shock. One 
can expect both a smeared shock and severe oscillations in the 
neighborhood of the shock. Hence, the grid work is constructed such 
that the ratio Af/A © is always near the maximum value for stabil- 
ity. The Mach numbers and angles of attack used in this study are 
such that no rapid expandions and, as a result, possible recompres- 
sion shocks are encountered in the 0 direction. Therefore, the 
increment A 0 is chosen such that the ratio a£/a 0 is as close 
to its maximum value as possible without having either too few or 
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an excessive number of grid points. The latter case leads to 
lengthy computation times. 

An estimate of the maximum mesh ratios and, as a result, the 
proper grid spacing can be made using the criteria given by Equations 
(43). Hence, the eigenvalues of the coefficient matricies of the 
gas dynamic equations of motion in the forebody coordinate system 
must be evaluated. The development of the expressions for the fore- 
body eigenvalues is presented in Appendix A. 

The gas dynamic equations of motion for the forebody are given 

i i 

by Equations (All) where the coefficient matricies A and B are 
defined in Equations (A12) - (A15). The five eigenvalues associated 
with each of the coefficient matricies are given by Equations (A16) 


- (A23) . 


The eigenvalues were determined numerically at all grid points 

in various preliminary solutions for a Mach number of 2 and angles 

of attack of 0 and 4 degrees. In all cases observed the triple 

repeated eigenvalues were the smallest. Thus, the largest of the 

eigenvalues given by Equations (A17), (A19), (A21) and (A23) 

represent the maximum for the corresponding matricies. For the 

» » 

wedge proper the maximum eigenvalues for the A and B matricies 
were near unity and for the tip half-cone 1.0 and 6.5 respectively. 


Since the eigenvalues of A for both the wedge proper and the 
cone tip have approximately the same maximum magnitude the corre- 
sponding A© increments for maximum stable values of the ratio 
Af/A© for both regions are very nearly the same. Therefore, 
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using a common A© in both regions is quite desirable from the 
stability viewpoint and certainly makes the computational process 
much simpler. 

On the other hand, in view of the maximum eigenvalues of the 
b' matricies, the A# increments for maximum stable values of the 
ratio A£/A0 should be near A® and 6.5 times A © for the wedge 
and cone tip regions respectively. 

The number of grid points in the 0 direction was set at 20 
with the A© increment chosen such that approximately 7 points are 
on the free stream side of the shock wave in the 2— D region of the 
flow field. The remaining grid points are in the shock layer region 
and sublayer. 

The A0 increments for the upper and lower wedge regions were 
chosen such that 10 grid points exist in each 9 = cons taint plane. 

At least 4 are outside the tip Mach cone and, as a result, in the 
2-D region of the flow field. The A0 increment for the cone tip 
region is chosen such that 13 grid points exist in each 9 = con- 
stant plane, 2 of which are common to the upper and lower wedge 
regions. 

The resulting computational plane grid system, similar to that 
depicted in Figure 23, has dimensions 20 by 31 which allowed 
reasonable solution convergence times. The numerical values 
associated with the A© and A0 are given in Table 1 for the 
two conditions investigated. 

It is noted upon examination of the numerical values for 
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Table 1. Computational plane grid spacing for the forebody flow 
region 



Upper 

wedge 

Cone 

tip 

Lower 

wedge 


CT= 0° 

„ „o 

a= 4 

a> o° 

cr= 4° 

a= o° 

CC= 4° 

Aj* 

(radians) 

0.1165 

0.1046 

0.2618 

0.2618 

0.1165 

0.1301 

A© 

(radians) 

0.04434 

0.04470 

0.04434 

0.04470 

0.04434 

0.04470 


the angular increments that the conditions for maximizing the mesh 
ratios are nearly satisfied. The worst violation occurs in the 
A0 increment for the upper and lower regions which are about 2.5 
times larger than the desired values. Decreasing the increments 
to their proper values, however, would require excessive grid 
points from a computation time viewpoint. Since no shock waves 
are encountered in the 0 direction, operation at the lower mesh 
ratios in these regions should not be prohibitive . 


Solution technique 

The integration is initiated in the 1 plane and proceeds 

20 steps in the £ direction at which place the solution is checked 
for convergence. If convergence has not been achieved the process 
is repeated until a solution is established. The increment is 
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chosen as large as possible without exceeding the stability bound 
which is determined experimentally by simply increasing until 
divergence occurs. 

In each £= constant plane (see Figure 23) the integration 
proceeds from the body surface to the free stream boundary and 
from the upper 2-D wedge boundary spanwise along the wing, around 
the cone tip and to the lower surface 2-D wedge boundary. 

Care must be taken in the integration process at the grid points 
in the wedge-cone interface planes since the equations of motion 
for the wedge proper and cone are not the same. Since the MacCormack 
predictor uses forward differences , the cone equations are used for 
the prediction step in the upper surface interface plane while the 
wedge equations are used in the lower surface interface plane. On 
the other hand, the opposite is true for the corrector step since 
the corrector utilizes backward differences. 

To insure that the predictor differences are always forward 
and the corrector differences backward the integration steps must 
always be in the positive directions. In view of the integration 
process described above, the positive directions for the y and z 
axes must be reversed in the lower surface wedge region. The only 
effect this has in the integration process is the reversal of signs 
associated with the v and w velocity components. 
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Forebodv numerical solutions 

One wing configuration is considered in this study with a fore 
body wedge half-angle of 7.5 degrees and, as a result, a 7.5 degree 
tip cone half-angle. The wing chord length is set at 2 units. 

The flow fields for a free stream Mach number of 2 are eval- 
uated at angles of attack of 0 and 4 degrees. The 0 degree case 
is used as a check to see if the numerical technique develops the 


proper flow field symmetry associated with the upper and lower 
surface flow regions. 

Flow field solutions for both the 0 and 4 degree angles of 
attack cases were obtained using the grid spacing defined in Table 
1. In both cases the near maximum integration step size was used. 

It was determined by means of a trial and error process. Succes- 
sively larger increments were tried until divergence occurred. 

The resulting A A increments for the 0 and 4 degree cases were 
0.04434 and 0.040 respectively. Table 2 shows the corresponding 

i 

mesh ratios A£/Aiz< and Af/A© used to obtain the solutions in the 
various regions of the forebody flow field. Also tabulated are the 
theoretical maximums based on the linear one-dimensional theory. 

It is noted that the ratio Aiff/A© is always within 80 per 
cent of the predicted maximum. This procedure should yield a 
reasonable solution in the shock vicinity. In the cross flow direc- 
tion the linear stability bound is exceeded in the cone region while 
on both the upper and lower wedge surfaces the ratio A^/A^ is 
between 30 and 40 per cent of the predicted maximum. As expected, 



Table 2. Theoretical and experimental mesh ratios used in the forebody solution 




a£/ai* 

/ M 


Af/A 

0 


Theoretical 

Value Per cent of 

Theoretical 

Value 

Per cent of 


maximum 

used 

theoretical 

maximum 

used 

theoretical 




maximum 



maximum 

0 

o 

II 

b 







Upper 

wedge 

1.017 

0.3805 

37.4 

1.039 

1.0 

96.25 

Tip 

cone 

0.1571 

0.1693 

107.8 

1.110 

1.0 

90.1 

Lower 

wedge 

1.017 

0.3805 

37.4 

1.039 

1.0 

96.25 

cr= 4 0 







Upper 

wedge 

0.9793 

0.3825 

39.1 

1.088 

0.895 

82.3 

Tip 

cone 

0. 1250 

0. 1526 

122.1 

1.111 

0.895 

80.6 

Lower 

wedge 

0.9430 

0.3075 

32.61 

1.057 

0.895 

84.7 
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increasing the increment causes instabilities in the crossflow 

direction in the cone region. 

The C£= 0° normal pressure distributions (normal to the wing 
chord plane) at various spanwise points in the mid-chord plane are 
shown in Figures 25 and 26 for the upper and lower wedge regions 
respectively. Figures 27 and 28 depict normal pressure distribu- 
tions in various meridian planes about the tip cone. The index 
parameter k defines the angular location of the 0 = constant planes 
in all regions. Figure 29 shows the relative location of the planes 
corresponding to the 31 k values. The index parameter j defines the 
angular location of the © s constant surfaces common to all regions. 
Table 3 lists the numerical values for © and 0 in all regions for 
both angles of attack considered. 

The distributions appear to be smooth with the bow shock ap- 
pearing in each distribution. It is quite well defined and usually 
contained in one to two intervals. As noted in Figure 30, the 
shock wave lies nearer the wing surface in the cone region. In the 
2-D wedge region the bow shock is at an angle of 36.71 degrees with 
respect to the wedge center line which decreases to 30.5 degrees in 
the cone region at the plane of symmetry (k = 16). The shock 
strength decreases from the wedge region to the cone region. 

The spanwise distribution of pressure along the body surface 
is shown in Figure 31. The body pressures decrease from the 2-D 
upper wedge value of 0.137 to a near constant value of about 0.1145 
around the cone surface then increases again to 0.137 at the lower 
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Figure 25. Forebody upper wedge region normal pressure distributions 
for the CT= 0° case 
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Figure 26. Forebody lower wedge region normal pressure distributions 
for the CC= 0° case 
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Forebody lower cone region normal pressure distributions 
for the CC = 0° case 
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Table 3. Numerical values for the angular orientations of the 
© = constant and 0 = constant planes defined by the 
indicies j and k respectively 


0 degrees © degrees 


k 

a= o° 

a= 4 ° 

J 

cr = o° 

a= 4° 

1 

-60.07 


1 

5.140 

4.939 

2 

- 53.40 

- 47.93 

2 

7.500 

7.500 

3 

- 46.72 

- 41.94 

3 

9.860 

10.06 

4 

- 40.05 

- 35.94 

4 

12.23 

12.62 

5 

- 33.37 

- 29.95 

5 

14.59 

15.18 

6 

- 26.70 

- 23.96 

6 

16.95 

17.74 

7 

- 20.02 

- 17.97 

7 

19.31 

20.30 

8 

- 13.35 

- 11.98 

8 

21.68 

22.86 

9 

- 6.670 

- 5.990 

9 

24.04 

25.42 

10 

0 

0 

10 

26.40 

27.99 

11 

15.0 

15.0 

11 

28.77 

30.54 

12 

30.0 

30.0 

12 - 

31.13 

33.11 

13 

45.0 

45.0 

13 

33.49 

35.67 

14 

60.0 

60.0 

14 

38.85 

38.23 

15 

75.0 

75.0 

15 

38.21 

40.79 

16 

90.0 

90.0 

16 

40.58 

43.35 

17 

105.0 

105.0 

17 

42.94 

45.91 

18 

120.0 

120.0 

18 

45.31 

48.47 

19 

135.0 

135.0 

19 

47.66 

51.03 

20 

150.0 

150.0 

20 

50.03 

53.60 

21 

165.0 

165.0 




22 

180.0 

180.0 




23 

- 6.67 

- 7.450 




24 

- 13.35 

- 14.91 




25 

- 20.02 

- 22.36 




26 

- 26.70 

- 29.81 




27 

- 33.37 

- 37.27 




28 

- 40.05 

- 44.72 




29 

- 46.72 

- 52.17 




30 

- 53.40 

- 59.63 




31 

- 60.07 

- 67.08 
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Figure 30. Forebody bow shock shape in the £ = 1.0 plane for the C C= 0° case 
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wedge 2-D region. Slight pressure oscillations occur along the 
body surface near the wedge-cone interface planes probably caused 
by the large values of the mesh ratio /\A//\0« These are the 
regions in which crossflow instabilities occur when the incre- 

ment is increased. 

In an attempt to demonstrate the reliability of the numerical 
technique used in this study, a comparison is made with flow fields 
developed by Babenko ^ 27 ^ for a circular cone with a half-angle of 
7.5 degrees at a Mach number of 2 at zero degrees angle of attack. 

The normal pressure distribution in the 90 degree meridian plane is 
plotted in Figures 27 and 28. The body surface distribution (which 
is constant) plotted around the cone surface from the 0 to the 180 
degree meridional planes is shown in Figure 31. The normal distri- 
butions in the 90 degree meridional plane as well as the body surface 
distribution for both the wedge— cone and the Babenko cone are similar 
in form with the numerical pressure values associated with the latter 
somewhat smaller. This is to be expected since the bow shocks asso- 
ciated with wedge shaped bodies are stronger them those associated 
with those of cones having the same vertex angle. 

In general, the numerical technique appears to generate the 
required flow field symmetry for the zero angle of attack case. 

The CC = 4° normal pressure distributions at various spanwise 
locations on the upper and lower wedge surfaces are shown in 
Figures 32 and 33 respectively. The distributions for the tip cone 
region in various meridional planes are depicted in Figures 34 and 35. 
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Figure 32. Forebody upper wedge region normal pressure distributions 
for the CC = 4° case 
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Figure 33. Forebody lower wedge region normal pressure distributions 
for the CT= 4° case 
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Figure 34. Forebody upper cone region normal pressure distributions 
for the CX = 4° case 
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Figure 35. Forebody lower cone region normal pressure distributions 
for the CC = 4° case 
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Again the distributions are very smooth with a well defined bow 
shock contained in one to two intervals. 

The normal pressure distribution in the 90 degree meridional 
plane (k = 16) for the Babenko cone at a 4 degree angle of attack 
is included in Figures 34 and 35. As expected, the pressured in 
the Babenko cone flow field are generally less than in the wedge- 
cone field. 

The shock shape in the £ = 1 plane for (X = 4° is shown in 
Figure 36. The shock angles with respect to the wedge center line 
in the upper and lower 2-D wedge regions have numerical values 
36.95 and 37.0 degrees respectively. In the cone region the bow 
shock lies generally near the body surface with the smallest center 
line shock angle (29.5 degrees) occurring in the 0 = 113 degree 
meridional plane on the lower surface. It is noted that the shock 
shapes for the CC - 0° and CC = 4° cases are quite similar. 

Although the center line shock angles in the upper and lower 
wedge regions are very nearly the same, the corresponding shock 
strengths are quite different as noted in the lateral body surface 
pressure distributions depicted in Figure 37. The surface pressure 
associated with the weaker upper surface shock is 0.1108 which de- 
creases only slightly around the cone then increases to 0.1681 in 
the lower 2-D wedge region containing the stronger bow shock. The 
slight pressure oscillations experienced in the CT = 0° case near 
the wedge-cone interface planes also appear in the (X = 4° data. 
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Region II Flow Field Analysis 

Afterbody geometry and coordinate systems 

That portion of the body contained between the mid-chord 
point and the trailing edge of the wing constitutes the afterbody 
which is a mirror image of the forebody. As with the forebody, 
the afterbody is separated into two parts. One part consists of 
the wedge proper and the other part the tip half-cone. 

The coordinate systems defining the afterbody tip cone and 
surrounding region is depicted in Figure 38. The Cartesian system 
is the same as that for the forebody. The system into which the 
equations of motion are cast consists of the ( A , 'Y , 0) system. 

The coordinate 0 denotes the angular orientation of the meridian 
plane containing the point ( £ , ~y , 0) and the axis of the tip half- 
cone. The coordinate 'Y (measured in the meridian plane) denotes 
the angle between the wing chord plane and a line in the meridian 
plane passing through the point ( £ , 'y , 0) and the circle given by 

y 2 + z 2 s (2 tan © c ) 2 (74) 

Note that the circle defined by Equation (74) is the intersection 
of the (y, z) plane and the extension of the tip cone. The 
coordinate £ represents simply the x position of the point 

. 7 . 0 ). 

The coordinate system used to define the wedge proper portion 
of the afterbody flow field is depicted in Figure 39. The coordi- 
nate 7 is the angle between the wing chord plane and the line that 
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contains the point (£, 'Y , 0) and the point y = 2 tan © w « The 
coordinate 0 represents the angle between the x axis and the plane 
containing the y axis and the point in question. Again, the 
coordinate ^ is simply the x coordinate. 

The grid system generated is defined in the tip region by 
the intersection of a set of 0 = constant meridian planes, a set 
of y a constant conical surfaces and a £ = constant plane. In 
the wedge proper region the grid system is defined by the inter- 
section of a set of 'Y = constant wedge planes, a set of 0 - 
constant planes and a £ = constant plane. It is noted that the 
afterbody surface is represented by a 'Y - constant surface, one 
reason for chosing the particular coordinate systems described. 

A typical grid network would be similar to that depicted in Figure 
23. Once again, the grid points in the wedge-cone interface 
planes are common to both the wedge and tip half-cone. 

i 

Afterbody equations of motion 

In order to perform the integration in the afterbody region 
the equations of motion given by Equations (37) must be transformed 
from the Cartesian form to the ( £.7. *> system. The trans- 
formation is presented in Appendix B. 

In the afterbody region the equations of motion for both the 
wedge proper and the tip half-cone are of the general form 


? 7 7 


0 


(75) 
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i i t 

The conservative variables E . F and G and the nonhomogeneous 


term H for the wedge proper are given by 


i 

E 


t 

F 


G* 


H 



-sin“yco^yE + cos 2 y F 
-sin^cos^E + cos 0G 


(76) 


(cos^y -sin'y - 2sin 2 0) E + 2sin0cosj2lG + 2sinycosyF 


and for the tip half-cone 


E = (£tany + 2tanO c ) E 


t 

F 


g' 


H 


2tan0 


= -sin 


inycosyE + cosy(cos)^F + sin^G) (tany+ — — — — ) 


-Fsin0 + Gcos0 

2tan6 _ (78) 

-(tany+ — -t — — ) (l-2cosy )E-2cosysin (cosjZfF+sinp^G) 


As in the forebody case, the variables E, F and G are the vector 
conservative variables associated with the Cartesian form of the 
equations of motion. 


Evaluation of the gas dynamic variables from conservative variables 
The same process is used to obtain the gas dynamic variables 
in the integration process over the afterbody region as is used in 
the forebody region. The first of Equations (76) is solved for the 
vector variable E" which is then incorporated in Equations (57) - (61) 
to evaluate v, w, u, jO and p. 
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Boundary conditions 

As in the forebody region, the afterbody boundary points must 
be specified in some fashion since the points along the edges of 
the grid are not integrated points. 

The grid spacing in the direction is chosen to be large 
enough that the outermost of the grid points are always in the 
free stream outside the shock layer. Hence, the gas dynamic 
variables along these boundaries are retained at the free stream 
values. 

Two methods are available for the specification of the variables 
along the 2-D afterbody wedge boundaries. Since the /\$ increments 
utilized in the afterbody region are the same as those used in the 
forebody integration, the boundary grid points as well as those in 
the three adjacent 0 = constant planes are all outside the tip 
Mach cone and in the 2-D flow region. Hence, the gas dynamic 
variables at the boundary may be set at the 2-D solution values that 
are known, a technique that was utilized in the forebody integration. 

A reflection technique may also be used at this boundary. The 
reflection occurs across the plane as shown in Figure 40. 


constant 

plane 


Figure 40. Geometry of the reflection technique 
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The gas dynamic variables at grid point 1 are set equal to those at 
grid point 3 as they must be if the points are in the 2-D region of 
the flow field. In effect, the technique generates its own 2-D 
solution. 

The reflection technique was adopted in this study. It pro- 
vides a means by which the integration technique can be tested. 

That is, if the solution in the 0 plane does not correspond to 
the known 2-D solution, the integration technique is not performing 
properly. 

The boundary conditions along the surface of the afterbody in 
each £ = constant plane are specified in the same fashion as on 
the forebody. A sublayer set of grid points is added and the 
modified Abbett technique is applied at the afterbody surface. It 
is noted here that the body entropies remain unaltered throughout 
the expansion at the mid-chord point of the wing. Hence, the 
entropy values are known from the forebody solution. 

Initial data 

The initial data required to start the afterbody integration 
are taken from the forebody solution. Since the grid points asso- 
ciated with the forebody and afterbody systems in a specified £ = 
constant plane do not coincide, point for point, some logical 
method must be used to transfer the data from the forebody to the 
afterbody grid system. In this study a scheme is utilized in 
which the four nearest forebody grid points for each afterbody 
grid point are located (see Figure 41). 
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0 - direction 



Figure 41. Data transfer technique from the forebody 
to the afterbody grid 

The data are then transferred by means of a simple linear inter- 
polation scheme in both the ff direction and *y direction. The 
specific £ = constant plane in which the data transfer occurs 
depends upon the particular case under consideration. Even though 
the data transfer occurs near the mid-chord the exact location of 
the initial data plane for the afterbody depends upon the step 
size used in the afterbody integration process. 

Choice of grid system 

The afterbody grid dimensions are established somewhat once 
the forebody system is defined in that the same 0 increments are 
used. This is to insure that the outer boundary points remain in 
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the 2-D region of the flow field. The increment £\Y is chosen 
such that at least four grid points^ always remain in the free stream 
region outside the shock layer throughout the afterbody region. 
Numerically, the values associated with AC/ for the afterbody and 
A© for the forebody are very similar in all cases. 


Solution technique 

The initial data plane for the afterbody is located such that 
the first integration step in the £ direction contains the mid- 
chord point at the half- interval . That is, the initial data plane 
in which the forebody data transfer occurs is located at £ = c/2 - 
where is the initial integration step size and c is the 

wing chord length. The initial integration step, then, is from the 
forebody region to the afterbody region. The integration then 
proceeds along the afterbody to the vicinity of the trailing edge. 
The magnitude of the integration step size is governed by the 

linear stability criteria. 

It was noted in earlier work that the gas dynamic equation 
coefficient matricies a' and B in the wedge proper region of the 


afterbody are identical in form to those associated with the fore- 
body. Hence, the same eigenvalues that dictate the maximum mesh 
ratios in the forebody wedge region now dictate the maximum of 
A£/A0 andA^/A X in the afterbody wedge proper region. These 
eigenvalues remain nearly constant over the afterbody region. 

For the tip cone, the eigenvalue structure is quite different 

i i 

in the afterbody region. The maximum eigenvalues for the A and B 
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matricies are given by Equations (A30) and (A32) respectively. 
Although the eigenvalues given by Equation (A30) are the same order 
of magnitude throughout the afterbody region, those defined by 
Equation (A32) become infinite as the trailing edge is approached. 
Hence, if constant increments based on initial data were used 

in the integration process one might expect a crossflow instability 
to develop since the maximum stable value for /\£/f\$ decreases 
rapidly near the trailing edge. In fact, the stable integration 
step size approaches zero near the trailing edge presenting dif- 
ficulties in the afterbody integration process. 

To insure that the crossflow instability does not occur the 
increment is set equal to the value dictated by the linear 

stability criteria given in Equations (43) at each step in the 
integration process. The value A^ is given by 


A0 

°^cone 


Xb' 


J 

max] 


(79) 


The integration moves quickly initially along the afterbody 
but slows down radically near the trailing edge. Furthermore, 
since all other eigenvalues remain nearly constant throughout the 
afterbody field, one can expect a lower quality flow field solu- 
tion near the trailing edge as a result of the very suboptimal mesh 
ratio operation, especially in the vicinity of the shock wave. 



108 


Afterbody numerical solutions 

Afterbody solutions are obtained using initial data from both 
the CC - 0° and CC = 4° forebody solutions at a Mach number of 2. 

The grid spacing in each £ = constant plane is quite similar 
to that of the forebody. The system contains 31 points in the 0 
direction and 20 points in the 'y direction. The £±0 increments 
for the afterbody region are the same as for the forebody region and 
are defined in Table 1. The A*y increments, common to all afterbody 
regions, are 3.0286 degrees and 3.0556 degrees for the GC = 0° and 
CC = 4° cases respectively. The angular orientation of the various 
0 = constant planes defined by the index parameter k are given in 
Table 3 and the angular orientation of the various *y - constant 
planes defined by the index parameter j are given in Table 4. 

An initial integration step size (^£) of 0.01 is used to 
move from the forebody region to the afterbody region. Hence, the 
initial data plane for the afterbody is at £ = 0.995 with the first 
integration advancing the data to the £ = 1.005 plane. The shoulder 
of the wing ( £ = 1.0) occurs at the half-interval of the first 
integration step. 

Subsequent integration step sizes are set to the maximum 
value dictated by the linear stability analysis given by Equation 
(79). The resulting numerical values for the mesh ratios (\^/ /\ 0 
and a£/aT are given in Tables 5 and 6. Included in the tables 
are the maximum values for the ratios in the various flow field 
regions as predicted by linear stability theory. Since the cone 


I 
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Table 4. Angular orientation of the various 'Y - constant planes 
in the afterbody grid system 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 
17 



*y (degrees) 



-10.53 

-7.50 

-4.47 

-1.44 

1.59 

4.62 

7.64 

10.67 

13.70 

16.73 

19.76 

22.39 

25.82 

28.84 

31.87 

34.90 

37.93 

40.96 

43.99 

47.02 


-10.56 

-7.50 

-4.44 

-1.39 

1.67 

4.72 

7.78 

10.83 

13.89 

16.95 

20.00 

23.06 

26.11 

29.17 

32.22 

35.28 

38.33 

41.39 

44.44 

47.50 


18 

19 

20 


Table 5. A comparison of the theoretical maximum values and the actual values of the 
mesh ratio a£/ao used in the afterbody integration 


Upper wedge 

( oc= 0°) 

Tip cone 

( a = o°) 

Lower wedge 

( cc = o°) 


Upper wedge 

( a = 4 °) 

Tip cone 

( a = 4°) 

Lower wedge 

( a = 4°) 


Theoretical 

maximum 


Value 

used 


Per cent of 
theoretical 
maximum 


£= 0.995 


£ = 0.995 

A = 1.97 

£ = 0.995 

• 

rl 

ii 

1.055 

2.185 

0.7523 

0.0205 

71.30 

0.940 

1.129 

2.208 

0.7523 

0.0205 

66.60 

0.93 

1.055 

2.185 

0.7523 

0.0205 

71.30 

1.16 

1.112 

1.997 

0.7815 

0.0397 

70.30 

1.85 

1.140 

2.011 

0.7815 

0.0397 

68.55 

1.97 

0.984 

2.183 

0.7815 

0.0397 

79.39 

1.82 
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Table 6. A comparison of the theoretical maximum values and the actual values of the 
mesh ratio A<^/A^ used in the afterbody integration 


Theoretical 

maximum 


Value 

used 


Per cent of 
theoretical 
maximum 


<£= 0.995 £= 1.97 £= 0.995 £ = 1.97 £ = 0.995 £ = 1 . 


Upper wedge 

( a = o°) 


1.015 1.933 0.3414 0.00932 33.64 0.48 


Tip cone 

( a = 0 °) 


0.1519 0.00415 0.1519 0.00415 100.0 100.0 


Lower wedge 

( a = o°) 


1.015 1.933 


0.3414 0.00932 33.64 


0.480 


Upper wedge 

( a = 4 °) 


1.098 1.904 


0.3986 0.02024 36.30 


1.060 


Tip cone 

( a = 4°) 


0.1592 0.00808 0.1592 0.00808 100.0 100.0 


Lower wedge 

( a= 4 °) 


1.053 2.174 


0.3204 0.01627 30.43 0.750 
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cross flow eigenvalues vary extensively throughout the integra- 
tion, the tables include values in the initial data plane ( £ = 
0 . 995 ) as well as the final integration plane ( £ = 1.97). 

As noted in Table 5, the ratio A£/Ay decreases from about 
70 per cent of its theoretical maximum initially to about 1 per 
cent near the trailing edge in order to satisfy the crossflow 
stability criteria. Hence, one can expect the flow field solution 
to degenerate somewhat in the vicinity of the shock as the trailing 
edge is approached. 

Normal pressure distributions at various spanwise locations on 
both the upper and lower wing surfaces at three different £ = con- 
stant locations along the afterbody chord are presented in Figures 
42 - 65. In particular. Figures 42 - 53 depict the normal pressure 
distributions for the C C = 0° angle of attack case in the £ = 1.2, 
1.5 and 1.76 planes along the afterbody. Figures 54 - 65 depict 
the normal pressure distributions for the d = 4° angle of attack 
in the £ = 1.25, 1.5 and 1.76 planes. The relative locations of 
the points on the wing surface (defined by the index parameter k) 
at which the normal distributions are displayed can be obtained in 
Figure 29 with the exact locations defined in Table 3. 

The exact 2-D solution for a symmetric double wedge afterbody 
is included on each set of upper and lower wedge surface pressure 
distributions. The exact solutions are used as the standard with 
which the numerical solutions in the 2-D region (k = 3) are com- 
pared to test the behavior of the numerical technique. 
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Figure 42. 
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Dimensionless Pressure 



Figure 43. Afterbody upper cone region normal pressure distribution in the 
£ = 1.2 plane for the CC = 0° case 
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Dimensionless Pressure 
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Dimensionless Pressure 
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Dimensionless Pressure 
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Figure 46. 


Afterbody upper wedge region normal p: 
in the £ = 1.5 plane for the CC = 0° 
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Dimensionless Pressure 
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Figure 
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Afterbody upper cone region normal pressure distributions 
in the fc = 1.5 plane for the CX = °° case 
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Dimensionless Pressure 



Figure 48. Afterbody lower wedge region normal pressure distributions 
in the t - 1.5 plane for the CT = 0 case 
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Dimensionless Pressure 
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Figure 49. Afterbody lower cone region normal pressure distributions 
in the C - 1.5 plane for the CC = 0° case 
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Dimensionless Pressure 
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Figure 51. 


Afterbody upper cone region normal pressure distributions 
in the £■ = 1.76 plane for the C£ = 0° case 
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Dimensionless Pressure 





Dimensionless Pressure 



Afterbody upper wedge region normal pressure distributions in 
the £ = 1.25 plane for the CT = 40 case 


Figure 54. 
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Dimensionless Pressure 
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Dimensionless Pressure 



'Y Direction Grid Index /V j 

Figure 56. Afterbody lower wedge region normal pressure distributions in 
the C = 1.25 plane for the GC = 4° case 
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Dimensionless Pressure 
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Figure 57. Afterbody lower cone region pressure distributions in the 
^ = 1.25 plane for the CC = 


4 case 


128 



Dimens io 



129 


Dimensionless Pressure 



Figure 59 . 


Afterbody upper cone region normal pressure distributions 
in the f = 1.5 plane for the CC = 4° case 
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Dimensionless Pressure 
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Dimensionless Pressure 
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Dimensionless Pressure 
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Dimensionless Pressure 
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Dimensionless Pressure 
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. Afterbody lower wedge region normal pressure distributions 
in the C = 1.76 plane for the CT = 4° case 


Figure 64 
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Dimensionless Pressure 
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The distributions, which now contain an expansion region as 
well as a bow shock, appear to be fairly smooth except in the 
vicinity of the shock where mild oscillations develop early in 
the integration process and become quite severe in the trailing 
edge region where the very low mesh ratios are used to satisfy 
the cross flow stability criteria. The shock waves in all dis- 
tributions, however, are quite well defined and usually contained 
in one to two intervals. The expansions are quite smooth and in 
the 2-D wedge regions follow the exact solutions quite closely. 

The numerical technique does, however, appear to overexpand the 
flow in both the wedge regions and cone regions on the body surface 
The 2-D wedge body pressure for the CT = 0° case is 0.0588, and for 
the CL = 4° case the upper and lower 2-D wedge body pressures are 
0.048 and 0.075 respectively. The numerical solution values are 
initially somewhat lower. However, as the integration proceeds to 
the trailing edge the 2-D body pressures increase again to nearly 
their proper values. The overexpansion also propagates into the 
flow field normal to the wing chord plane. 

With the exception of the overexpansion situation near the 
body surface, the numerical technique generates a 2-D solution that 
compares favorably with the exact solution. In addition, the pres- 
sure distributions for the CT = 0° case indicate that the proper 
flow field symmetry is generated. 
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Region III Flow Field Analysis 

Afterb ody geometr y, c oordinate system and equations of motion 

That portion of the flow field aft of the trailing edge of 
the wing is termed the wake. This region, void of any solid 
boundaries, is bounded on three sides by the free stream and on 
the fourth side by the 2-D flow field associated with the wake of 
a symmetric double Wedge. 

The coordinate system used to define the wake flow field is 
a Cartesian system with the origin centered at the forebody cone 
apex. The same frame was used to describe the body geometry 
associated with the forebody and the afterbody. The grid in each 
x = constant plane is defined by the intersection of a set of y = 
constant and z = constant planes, resulting in a rectangular 
arrangement of mesh points. 

The equations of motion in the Cartesian frame are given by 
Equations (38) with the conservative variables IT, F and G given in 
Equations (39) - (41). 

Initial and boundary conditions 

The boundary conditions applied at the outer extremes of the 
grid are similar to those used in the afterbody integration. The 
grid dimensions are chosen large enough such that the outer grid 
points always lie in the free stream along the upper and lower 
grid boundaries as well as on the boundary off the wing tip. Along 
the 2-D flow region boundary the reflection technique is used to 
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determine the flow variables. Hence, as long as the flow along 
this boundary remains two-dimensional the numerical solution may 
be compared with the known 2-D exact solution in an attempt to 
evaluate the performance of the numerical technique. It is noted 
here that if the tip Mach cone intersects this boundary the flow 
will no longer be two-dimensional and the 2-D exact solution no 
longer applies. However, the 3-D solution generated by the numer- 
ical technique in this case is still legitimate with the z = constant 
plane about which the reflection occurs representing the center of 
the 3-D wing. 

The initial data used in the wake integration are generated 
in the afterbody integration and are contained in the final integra- 
tion plane of the afterbody. Since the grid points associated with 
the afterbody and wake systems do not, in general, coincide the 
data must be transferred from the afterbody grid to the wake grid. 

The linear interpolation scheme used in the forebody-afterbody 
data shift is used here. 

Grid system and solution technique 

The wake grid in each x = constant plane contains 40 mesh 
points in both the y and z directions. The grid is positioned in 
the wake region with the j = 20 and 21 mesh points at y = /\y/2 
and - Ay/2 respectively (see Figure 66). Hence, half the grid 
points lie above the wing chord plane and half below. 
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Figure 66. A portion of the wake grid system showing the 
relative location of the wing 


The numerical values for /\y and /\z used in the wake inte- 
gration are 0.195 and 0.176 respectively for the CC = 0° case and 
0.195 and 0.1584 respectively for the Cf = 4° case. The incremental 
magnitudes were chosen such that in the wake initial data plane the 
3-D portion of the wing flow field is contained within a region 
bounded on the lower and upper sides by grid points with indicies 
j = 10 and 30 respectively and on the left and right by grid points 
with indicies k = 10 and 30 respectively^ Hence, the outermost 10 
grid points in the system lie either in the free stream or in the 
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2-D double-wedge flow field. However, as the integration proceeds 
into the wake the 3— D flow field propagates throughout the grid to 
the grid boundaries. 

The afterbody integration was terminated at x = 1.97. Hence, 
the flow variables at the grid points in the x = 1.97 plane provide 
the initial data for the wake integration. The data are transferred 
from the afterbody grid to the wake grid, and a single integration 
step of magnitude Ax — 0.06 generates the data in the first x — 
constant plane in the wake behind the wing. It is noted that the 
wing trailing edge is located at the half— interval of the first 
integration step. 

Subsequent integration step sizes were chosen to be the maximum 
allowed as predicted by the linear stability theory. The eigenvalues 
associated with the equations of motion for the wake are developed 
in Appendix A. Equations (A6) and (A9) represent the maximum eigen- 
values of the A and B matricies respectively and, as a result, are 
used in conjunction with Equations (43) to determine the maximum 
step size Ax. For both cases considered the integration was 
terminated well before the shock wave intersects the grid boundaries 
which occurs at about one chord length behind the wing. 

Numerical solutions 

Figures 67 and 68 show various pressure distributions normal 
to the wing chord plarie for the CC = 0° case in the x = 2.57 and 
x = 3.0 planes respectively. Figures 69 and 70 show similar 



Dimensionless Pressure 



Figure 67. Pressure distributions normal to the wing chord plane for the 
CT = 0° case in the x = 2.57 plane 
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Dimensionless Pressure 



. Pressure distributions normal to the wing chord plane for the 
CC = 0° case in the x = 3.0 plane 


Figure 68 
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Dimensionless Pressure 
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Figure 69. Pressure distributions normal to the wing chord plane for the 
CC = 4° case in the x = 2.12 plane 
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Dimensionless Pressure 
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Figure 70. Pressure distributions normal to the wing chord plane for the 
C C ~ 4° case in the x = 2.6 plane 
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distributions for the CC ~ 4 case in the x = 2.12 and x = 2.6 planes 
respectively. Each contains an exact 2-D solution for comparison. 

It is noted again here that the index values k = 1 , 20 and 40 define 
grid points in the 2-D double-wedge wake flow region, the wing tip 
and free stream region off the wing tip respectively. 

In moving along the pressure distributions point to point from 
j = 1 to j = 40 one encounters a bow shock, an expansion and a re- 
compression shock all associated with the lower surface flow, then 
another recompression shock, an expansion and a bow shock associated 
with the upper surface flow. 

The distributions appear to be smooth with well defined 
shocks as in the forebody and afterbody solutions. The numerical 
2-D solution (k = 1) appears to agree quite well with the exact 
2-D solution. The bow shocks have proper strengths as well as loca- 
tions. However, the strengths of the recompression shocks as pre- 
dicted by the numerical method appear to be excessive. For the 
a = o° case the pressure behind the recompression shock should 
be 0.0935 with the worst numerical solution yielding 0.101, re- 
sulting in an 8.02 per cent error. For the CC = 4° case the exact 
and worst numerical values are 0.94 and 0.1075 respectively, 
yielding a 14.4 per cent error. The recompression shock locations 
in the 2-D region seem to be well predicted by the numerical method 


in all cases. 
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RECOMMENDATIONS FOR FURTHER STUDY 

One of the major difficulties associated with the wake region 
analysis is that the Cartesian coordinate system used to generate 
the wake grid system severely limits the distance into the wake 
that the integration can proceed. Once the shock intersects the 
grid boundary the integration must be terminated for lack of 
proper boundary conditions. To proceed farther into the wake the 
grid system must be enlarged. This increases the computer storage 
requirements and quickly becomes prohibitive in terms of computer 
capacity. 

This problem can be eliminated by changing the wake coordinate 
system to a conical system with an origin at the apex of the fore- 
body cone. Then a grid system can then be developed such that the 
3-D flow field is contained within the grid boundaries regardless 
of distance behind the wing. 

The analysis could also be extended to include a variety of 
wing cross-sections such as parabolic or circular arc sections, 
cambered as well as uncambered. In addition, a wider variety of 
tip geometries could be investigated. 

Although exact solutions are available for comparison with 
the numerical solution in the 2-D regions of the flow field no 
such standards exist for the 3-D flow regions. Hence, the per- 
formance of the numerical technique can only be extrapolated from 
the 2-D results. It would be quite desirable to have experimental 
evidence to verify the numerical results in the 3-D regions. 
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Although in so far as the author knows no experimental investiga- 
tions have been made on the particular wing configuration used in 
this analysis, experimental studies have been made on rectangular 

/ OQ \ 

wings of other cross sections. In particular, Davis conducted 

experimental tests on various nonlifting rectangular planform wings 
with parabolic cross-sections. The numerical technique used in this 
investigation could be easily adapted to the Davis wing configura- 
tions. Hence, the 3-D region performance of the numerical technique 

could be checked. 

The recommendations suggested above are based on the experience 
gained in applying the MacCormack technique to the wing configura- 
tion used in this study. The method yielded satisfactory results 
in most cases and, as a result, appears to represent a powerful 
tool which could be used to investigate flow fields about various 
other configurations. 
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APPENDIX A 

Eigenvalue Evaluation for Cartesian System 

The stable range of mesh ratios contained in the finite 
difference equations can be theoretically. predicted using am- 
plification matrix theory. Application of this theory requires an 
evaluation of the eigenvalues associated with the coefficient 
matricies of the gas dynamic equations of motion. 

Kutler^ 5 ^ shows that the steady equations in a Cartesian 
frame can be written in the form 


U +AU +BU +C=0 
X y z 


where U represents a vector of state variables 



P 


[P) 


(Al) 


(A2) 


and where A and B are square matrices given by 
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(A4) 


In the development of the matrices A and B the energy equation 
was used in the differential form 

= 0 (A5) 


q * 


Vp - ( 2p/a/0) s V<> 
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For the A matrix, three of the eigenvalues are identical and 
g iven by 


K. 2, 3 * v/u 


while the remaining two are given by 


(A6) 


\ -uv 1 cVu 2 + v 2 - c 2 

A 4, 5 


2 2 
u - c 


(A7) 


In a similar fashion, three of the eigenvalues of the B matrix 
are identical and given by 


X l, 2, 3 = " /U 


whereas the remaining two are given by 


(A8) 


+ -v / 2 2 2 

X - ~ uw - C + w - c 

4, 5 2 2 

* u - c 


(A9) 


Forebody Eigenvalues 


The steady gas dynamic equations for the forebody wedge proper 
and the tip half-cone can be developed by a simple transformation 
of independent variables in Equations (Al) using the Jacobian 
elements in Table B2. In terms of the ( 9, f6) system the 

equations become 



(A10) 


+ A <U. £ y ♦ U e S y ♦ ll/ y ) 

+b ( u t L ♦ Vz - y j * c * ° 

Substituting for the Jacobian elements yields a set of equations 
given by 


+ A U Q + B U0 + C = 0 
where, for the wedge proper 


. » A cos 9 I sin 9 Cos 9 
A = 


? 


„ ' B cos^ <j/> I sin $ cos $ 

€ ~ f 

and for the tip half-cone 


2 2 

’ _ Acos 9cos0 Bcos 9sin0 Isin9cos9 
_ Bcos#cos9 Asinj#cos9 




(All) 


(A12) 

(A13) 


(A14) 

(A15) 


The matrix I is the identity matrix. In their fully expanded form 
showing all the elements the matrices defined by Equations (A12) - 
(A15) are presented on the following four pages, Equations (A12a) - 


(A15a) . 

Each of the four matrices contain at least one row or column 
in which all but one of the elements are zero. Hence, one of the 
eigenvalues in each matrix is immediately available. In all cases 
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it can be shown that the first eigenvalue is a triple repeated 
eigenvalue. The remaining two eigenvalues are then determined 
by solving the appropriate quadratic expression. 

I 

For the wedge proper the eigenvalues of the A matrix are 
given by Equations (A16) - (A17). 
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1, 


2 , 



-sinQcos© vcos 9 
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2 2 
V c 
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For the B matrix the eigenvalues are 


X -singfcosgf wcos 

1, 2, 3 = £ + £ u 

X - ~ s i n 0 cos 0 uwcos 0 
4 - 5 = £ ' £<c 2 -u 2 > 


c cos^ 


£ (c -u ) L 
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u - 


given by 



» 

For the tip half-cone the eigenvalues of the A matrix are 
given by 


2 2 

X -sinQcos© vcos Qcosgf wcos 9s in# 

1, 2, 3 = £ 




4, 5 


2 2 2 
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(A16) 

(A17) 

(A18) 

(A19) 

(A20) 


(A21) 
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and for the B matrix 


^ _ -v cos 9 sin 0 + w cos 8 cos 0 

1, 2, 3 = U 


X — uw cosQ cos# 
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Afterbody Eigenvalues 

For the afterbody, the transformed equations of motion using 
the Jacobian elements of Table B4 in the Cartesian equations of 
motion, Equations (Al), are given by 


U + A U + B Uj + C = 0 

c 7 

Where for the wedge proper 

. 1 A cos^V I sinV cos"/ 
A 

* _ B cos $ I sin 0 cos 0 

' * ' £ 

and for the tip half -cone 


1 A 
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(A24) 


(A25) 

(A26) 


(A27) 


(A28) 
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The A* and B* matricies for the afterbody wedge proper are the 
same as those for the forebody with the *y parameter replacing the 

f I 

Q parameter. Hence, the matricies A and B for the afterbody are 
given by Equations (A12a) and (A13a) respectively and the resulting 
eigenvalues are given by Equations (A16) - (Ai9) with the *y param- 
eter replacing the 9 parameter. 

1 * . 

The fully expanded form of the A and B matricies for the 
tip half-cone, as defined in Equations (A27) and (A28), are given 
by Equations (A27a) and (A28a) on the following two pages. The 
matricies are similar to the corresponding matricies in the forebody 
region in that three of the A eigenvalues for the tip cone, (A29), 
are triple repeated. The remaining two eigenvalues are given by 
Equation (A30) . For the b' matrix the eigenvalues are given by 
Equations (A31) and (A32) . 
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APPENDIX B - FLOW EQUATION TRANSFORMATIONS 
Introduction 

Integration of the equations of motion via MacCormack's tech- 
nique requires that the equations be cast in a conservative form 
using the independent variables associated with the various coordinate 
systems describing the body and' flow field. Hence, the Cartesian 
conservative form presented in Equations (37) must be transformed 
in Regions I and II to the various systems defining the forebody 
and afterbody. The Jacobian elements of the transformations as well 
as the transformed equations of motion and resulting conservative 
forms are presented in subsequent sections. 

Forebody Equations of Motion 

The coordinate systems used to describe the forebody are given 
in Figures 21 and 22. Although the coordinate labels ( £ > *’ e) 
are identical their definitions are, in general, different. 

In the Cartesian system it has been shown that the equations 
of motion can be written in the form 

E + F + G = 0 (Bl) 

x y z 

where the vector conservative variables E, F, and G are given 
by Equations (39), (40) and (41). Noting that 



167 


x = x(£, 0, ©) 

y = y( £ » 0 > ©) ( B2 ) 

2 = z(£ > 0 , ©) 

the equations of motion (Al) can be written as 


C + E„/ 0 + E. 0 

® x 0 'x © > 


+ F. £ + F . 0 + F © (B3) 

£ ^y 0 y © y 

+ g. C + g . 0 + g? e = o 

t 0 Z © z 

where the terms C > C > £ > 0 > 0 t 0 > & > ® > an d ® represent 

3 x 3 y 3 z x y z x y z 

the Jacobian elements of the transformation. 

The transformation equations for the wedge proper and the half- 
cone are given in Table Bl. The resulting Jacobian elements are 
given in Table B2. 

Substitution of the Jacobian elements of Table B2 into 
Equations (Al) yields for the wedge proper 


[£e 1 + T -sin©cos©E + cos 2 0F 1 

* L 

♦ [-sinjfcos# + cos' '0C J 
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0 


(B4) 


[ 2 2 2 — M 1 
(-sin ©-sin 0+cos 0)B + 2sin0cos©F + 2sin0cos^G = 0 


and for the half -cone 
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^£tan©Ej + J”-Esin 2 © + Fsin©cos0 + Gsin©cos©sin0 J 
+ -Fsin0 + Gcos0 I 

L J0 


“ 2 — 2 

+ (2sinQcos© - tan©)E + 2Fcos0sin 9 + 2Gsin 9sin0 


Equations (B4) and (B5) represent the conservative form of the 
equations of motion in the forebody region. 


Table Bl. Forebody transformation equations 


Wedge proper 


Half -cone 


tan 1 z/x 


tan 1 z/y 


tan” 1 y/x 


tan x 
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transformation 


Half-cone 
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Afterbody Equations of Motion 

The coordinate systems used to describe the afterbody are 
given in Figures 38 and 39. The equations of motion in the 
Cartesian form are given by Equation (Bl) and in the ■(£ , 0 , 'Y ) 
system by Equation (B3) with the 9 variable replaced by the 
variable. The transformation equations 

x = x( 0,y ) 

y = y(£ » 0,y ) ( B6 ) 

z = ■(£. 0,y) 

are given in Table B3 for the wedge proper and the tip half-cone 

with the resulting Jacobian elements £ . f L, . 0„. 0,. 

x y z x y z 

y , y , and y given in Table B4. 
x r y 'z 


Table B3. Coordinate transformation equations for the afterbody 
region equations of motion 



Wedge proper 

Half -cone 


X 
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tan -1 (z/x) 

tan ^(z/y) 
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2 2 
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“ lyj^-z2 J 

L x J 
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Table B4. Jacobian elements for the afterbody transformation 


Wedge proper 


Half -cone 


Substitution of the Jacobian elements of Table B4 into Equa- 
tions (Al) yields for the afterbody wedge proper 

•sin0cosJ2^E + cos 


<f E + -si] 


-sin cy cosy 


■0 


E + cos^y f| + j^-2si 

r 
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and for the tip half-cone 
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